<p><b>dwj07@fsu.edu</b> 2013-01-31 10:05:58 -0700 (Thu, 31 Jan 2013)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
First cut of OpenMP on element loops.<br>
This includes the bug fixes related to scratch variables, and the ability to create scratch super arrays.<br>
<br>
Equation of state is not OpenMP'd but the rest of the loops are bit reproducible comparing OpenMP to flat MPI.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/openmp_elements/Makefile
===================================================================
--- branches/ocean_projects/openmp_elements/Makefile 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/Makefile 2013-01-31 17:05:58 UTC (rev 2403)
@@ -131,7 +131,8 @@
"DEBUG = $(DEBUG)" \
"SERIAL = $(SERIAL)" \
"USE_PAPI = $(USE_PAPI)" \
- "CPPFLAGS = $(MODEL_FORMULATION) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" )
+ "CPPFLAGS = $(MODEL_FORMULATION) -DUNDERSCORE -m64 $(FILE_OFFSET) $(ZOLTAN_DEFINE)" \
+ "OPENMP_FLAG = -fopenmp" )
g95:
( $(MAKE) all \
@@ -277,6 +278,15 @@
TAU_MESSAGE="TAU Hooks are off."
endif
+ifeq "$(OPENMP)" "true"
+ FFLAGS+=$(OPENMP_FLAG)
+ CFLAGS+=$(OPENMP_FLAG)
+ LDFLAGS+=$(OPENMP_FLAG)
+ OPENMP_MESSAGE="OpenMP is on."
+else
+ OPENMP_MESSAGE="OpenMP is off."
+endif
+
ifneq ($(wildcard $(NETCDF)/lib/libnetcdff.*), ) # CHECK FOR NETCDF4
LIBS += -lnetcdff
endif # CHECK FOR NETCDF4
@@ -316,6 +326,7 @@
@echo ""
@echo $(DEBUG_MESSAGE)
@echo $(SERIAL_MESSAGE)
+ @echo $(OPENMP_MESSAGE)
@echo $(PAPI_MESSAGE)
@echo $(TAU_MESSAGE)
clean:
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/Registry 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/Registry 2013-01-31 17:05:58 UTC (rev 2403)
@@ -323,3 +323,27 @@
var persistent integer edgeSignOnCell ( maxEdges nCells ) 0 - edgeSignOnCell mesh - -
var persistent integer edgeSignOnVertex ( maxEdges nVertices ) 0 - edgeSignOnVertex mesh - -
var persistent integer kiteIndexOnCell ( maxEdges nCells ) 0 - kiteIndexOnCell mesh - -
+
+% Scratch fields, for openmp
+var scratch real uTempEdges ( nEdges ) 0 - uTempEdges scratch - -
+var scratch real delsq_u ( nVertLevels nEdges ) 0 - delsq_u scratch - -
+var scratch real delsq_vorticity ( nVertLevels nVertices ) 0 - delsq_vorticity scratch - -
+var scratch real delsq_divergence ( nVertLevels nCells ) 0 - delsq_divergence scratch - -
+var scratch real uh ( nVertLevels nEdges ) 0 - uh scratch - -
+var scratch real high_order_horiz_flux ( nVertLevels nEdges ) 0 - high_order_horiz_flux scratch - -
+var scratch real tracer_new ( nVertLevels nCells ) 0 - tracer_new scratch - -
+var scratch real tracer_cur ( nVertLevels nCells ) 0 - tracer_cur scratch - -
+var scratch real upwind_tendency ( nVertLevels nCells ) 0 - upwind_tendency scratch - -
+var scratch real inv_h_new ( nVertLevels nCells ) 0 - inv_h_new scratch - -
+var scratch real tracer_max ( nVertLevels nCells ) 0 - tracer_max scratch - -
+var scratch real tracer_min ( nVertLevels nCells ) 0 - tracer_min scratch - -
+var scratch real flux_incoming ( nVertLevels nCells ) 0 - flux_incoming scratch - -
+var scratch real flux_outgoing ( nVertLevels nCells ) 0 - flux_outgoing scratch - -
+var scratch real high_order_vert_flux ( nVertLevelsP1 nCells ) 0 - high_order_vert_flux scratch - -
+var scratch real delsq_temperature ( nVertLevels nCells ) 0 - delsq_temperature scratch delsq_tracers dynamics
+var scratch real delsq_salinity ( nVertLevels nCells ) 0 - delsq_salinity scratch delsq_tracers dynamics
+var scratch real delsq_tracer1 ( nVertLevels nCells ) 0 - delsq_tracer1 scratch delsq_tracers testing
+var scratch real drhoTopOfCell ( nVertLevelsP1 nCells ) 0 - drhoTopOfCell scratch - -
+var scratch real drhoTopOfEdge ( nVertLevelsP1 nEdges ) 0 - drhoTopOfEdge scratch - -
+var scratch real du2TopOfCell ( nVertLevelsP1 nCells ) 0 - du2TopOfCell scratch - -
+var scratch real du2TopOfEdge ( nVertLevelsP1 nEdges ) 0 - du2TopOfEdge scratch - -
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_equation_of_state_jm.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_equation_of_state_jm.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_equation_of_state_jm.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -264,6 +264,7 @@
enddo
endif
+ !!$OMP DO PRIVATE(k, SQ, TQ, SQR, T2, WORK1, WORK2, RHO_S, WORK3, WORK4, BULK_MOD, DENOMK)
do iCell=1,nCells
do k=1,maxLevelCell(iCell)
SQ = max(min(tracers(indexS,k,iCell),smax),smin)
@@ -308,6 +309,7 @@
end do
end do
+ !!$OMP END DO
deallocate(pRefEOS,p,p2)
end subroutine ocn_equation_of_state_jm_rho!}}}
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_gm.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_gm.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_gm.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -52,16 +52,20 @@
if (config_vert_grid_type .EQ. 'isopycnal') then
+ !$OMP DO PRIVATE(k)
do iEdge = 1, nEdges
do k = 1, maxLevelEdgeTop(iEdge)
uBolusGM(k,iEdge) = hEddyFlux(k,iEdge)/h_edge(k,iEdge)
end do
end do
+ !$OMP END DO
else
! Nothing for now for all other grid types (zlevel, zstar, ztilde)
+ !$OMP WORKSHARE
uBolusGM(:,:) = 0.0
+ !$OMP END WORKSHARE
end if
@@ -88,9 +92,12 @@
nEdges = grid % nEdges
+ !$OMP WORKSHARE
hEddyFlux(:,:) = 0.0
+ !$OMP END WORKSHARE
if (config_vert_grid_type .EQ. 'isopycnal') then
+ !$OMP DO PRIVATE(cell1, cell2, k)
do iEdge = 1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -98,6 +105,7 @@
hEddyFlux(k,iEdge) = -config_h_kappa * (h(k,cell2) - h(k,cell1)) / dcEdge(iEdge)
end do
end do
+ !$OMP END DO
else
!Nothing for now for all other grid types (zlevel, zstar, ztilde)
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_mpas_core.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_mpas_core.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -344,27 +344,33 @@
integer :: ierr
! Eventually, dt should be domain specific
+ !$OMP PARALLEL
+
+ !$OMP BARRIER
+ !$OMP MASTER
dt = config_dt
-
currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
call mpas_get_time(curr_time=currTime, dateTimeString=timeStamp, ierr=ierr)
+
write(0,*) 'Initial time ', trim(timeStamp)
-
if (config_write_output_on_startup) then
call ocn_write_output_frame(output_obj, output_frame, domain)
endif
+
block_ptr => domain % blocklist
-
do while(associated(block_ptr))
call ocn_time_average_init(block_ptr % state % time_levs(1) % state)
block_ptr => block_ptr % next
end do
+ !$OMP END MASTER
+ !$OMP BARRIER
! During integration, time level 1 stores the model state at the beginning of the
! time step, and time level 2 stores the state advanced dt in time by timestep(...)
itimestep = 0
do while (.not. mpas_is_clock_stop_time(clock))
-
+ !$OMP BARRIER
+ !$OMP MASTER
itimestep = itimestep + 1
call mpas_advance_clock(clock)
@@ -379,7 +385,11 @@
end do
call mpas_timer_start("time integration", .false., timeIntTimer)
+ !$OMP END MASTER
+ !$OMP BARRIER
call mpas_timestep(domain, itimestep, dt, timeStamp)
+ !$OMP BARRIER
+ !$OMP MASTER
call mpas_timer_stop("time integration", timeIntTimer)
! Move time level 2 fields back into time level 1 for next time step
@@ -388,7 +398,10 @@
call mpas_shift_time_levels_state(block_ptr % state)
block_ptr => block_ptr % next
end do
-
+ !$OMP END MASTER
+
+ !$OMP BARRIER
+ !$OMP MASTER
if (mpas_is_alarm_ringing(clock, outputAlarmID, ierr=ierr)) then
call mpas_reset_clock_alarm(clock, outputAlarmID, ierr=ierr)
! output_frame will always be > 1 here unless it was reset after the
@@ -421,8 +434,10 @@
call mpas_output_state_for_domain(restart_obj, domain, 1)
call mpas_output_state_finalize(restart_obj, domain % dminfo)
end if
-
+ !$OMP END MASTER
+ !$OMP BARRIER
end do
+ !$OMP END PARALLEL
end subroutine mpas_core_run!}}}
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_restoring.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_restoring.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_restoring.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -126,6 +126,7 @@
invSalinity = 1.0 / (salinityTimeScale * 86400.0)
k = 1 ! restoring only in top layer
+ !$OMP DO
do iCell=1,nCellsSolve
tend(indexT, k, iCell) = tend(indexT, k, iCell) - h(k,iCell)*(tracers(indexT, k, iCell) - temperatureRestore(iCell)) * invTemp
tend(indexS, k, iCell) = tend(indexS, k, iCell) - h(k,iCell)*(tracers(indexS, k, iCell) - salinityRestore(iCell)) * invSalinity
@@ -136,6 +137,7 @@
! / (config_restoreT_timescale * 86400.0)
enddo
+ !$OMP END DO
!--------------------------------------------------------------------
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tendency.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tendency.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -102,18 +102,21 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tend_h(tend, s, grid)!{{{
+ subroutine ocn_tend_h(tend, s, grid, scratch)!{{{
implicit none
type (tend_type), intent(inout) :: tend !< Input/Output: Tendency structure
type (state_type), intent(in) :: s !< Input: State information
type (mesh_type), intent(in) :: grid !< Input: Grid information
+ type (scratch_type), intent(inout) :: scratch !< Input/Output: Scratch structure
real (kind=RKIND), dimension(:,:), pointer :: h_edge, wTop, tend_h, uTransport
integer :: err
+ !$OMP SINGLE
call mpas_timer_start("ocn_tend_h")
+ !$OMP END SINGLE
uTransport => s % uTransport % array
wTop => s % wTop % array
@@ -124,7 +127,10 @@
!
! height tendency: start accumulating tendency terms
!
+
+ !$OMP WORKSHARE
tend_h = 0.0
+ !$OMP END WORKSHARE
!
! height tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( hu)
@@ -134,18 +140,24 @@
!
! QC Comment (3/15/12): need to make sure that uTranport is the right
! transport velocity here.
+ !$OMP SINGLE
call mpas_timer_start("hadv", .false., thickHadvTimer)
+ !$OMP END SINGLE
call ocn_thick_hadv_tend(grid, uTransport, h_edge, tend_h, err)
+ !$OMP SINGLE
call mpas_timer_stop("hadv", thickHadvTimer)
!
! height tendency: vertical advection term -d/dz(hw)
!
call mpas_timer_start("vadv", .false., thickVadvTimer)
+ !$OMP END SINGLE
call ocn_thick_vadv_tend(grid, wtop, tend_h, err)
+ !$OMP SINGLE
call mpas_timer_stop("vadv", thickVadvTimer)
call mpas_timer_stop("ocn_tend_h")
+ !$OMP END SINGLE
end subroutine ocn_tend_h!}}}
@@ -162,13 +174,14 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tend_u(tend, s, d, grid)!{{{
+ subroutine ocn_tend_u(tend, s, d, grid, scratch)!{{{
implicit none
type (tend_type), intent(inout) :: tend !< Input/Output: Tendency structure
type (state_type), intent(in) :: s !< Input: State information
type (diagnostics_type), intent(in) :: d !< Input: Diagnostic information
type (mesh_type), intent(in) :: grid !< Input: Grid information
+ type (scratch_type), intent(inout) :: scratch
real (kind=RKIND), dimension(:,:), pointer :: &
h_edge, h, u, rho, zMid, pressure, &
@@ -179,7 +192,9 @@
integer :: err
+ !$OMP SINGLE
call mpas_timer_start("ocn_tend_u")
+ !$OMP END SINGLE
u => s % u % array
rho => s % rho % array
@@ -204,32 +219,42 @@
! velocity tendency: start accumulating tendency terms
!
! mrp 110516 efficiency: could remove next line and have first tend_u operation not be additive
+ !$OMP WORKSHARE
tend_u(:,:) = 0.0
+ !$OMP END WORKSHARE
!
! velocity tendency: nonlinear Coriolis term and grad of kinetic energy
!
+ !$OMP SINGLE
call mpas_timer_start("coriolis", .false., velCorTimer)
+ !$OMP END SINGLE
+
call ocn_vel_coriolis_tend(grid, Vor_edge, h_edge, u, ke, tend_u, err)
+ !$OMP SINGLE
call mpas_timer_stop("coriolis", velCorTimer)
!
! velocity tendency: vertical advection term -w du/dz
!
call mpas_timer_start("vadv", .false., velVadvTimer)
+ !$OMP END SINGLE
call ocn_vel_vadv_tend(grid, u, h_edge, wtop, tend_u, err)
+ !$OMP SINGLE
call mpas_timer_stop("vadv", velVadvTimer)
!
! velocity tendency: pressure gradient
!
call mpas_timer_start("pressure grad", .false., velPgradTimer)
+ !$OMP END SINGLE
if (config_pressure_type.eq.'MontgomeryPotential') then
call ocn_vel_pressure_grad_tend(grid, MontPot, zMid, rho, tend_u, err)
else
call ocn_vel_pressure_grad_tend(grid, pressure, zMid, rho, tend_u, err)
end if
+ !$OMP SINGLE
call mpas_timer_stop("pressure grad", velPgradTimer)
!
@@ -238,7 +263,9 @@
! strictly only valid for config_h_mom_eddy_visc2 == constant
!
call mpas_timer_start("hmix", .false., velHmixTimer)
- call ocn_vel_hmix_tend(grid, divergence, vorticity, viscosity, tend_u, err)
+ !$OMP END SINGLE
+ call ocn_vel_hmix_tend(grid, divergence, vorticity, viscosity, tend_u, scratch, err)
+ !$OMP SINGLE
call mpas_timer_stop("hmix", velHmixTimer)
!
@@ -248,18 +275,25 @@
! know the bottom edge with nonzero velocity and place the drag there.
call mpas_timer_start("forcings", .false., velForceTimer)
+ !$OMP END SINGLE
call ocn_vel_forcing_tend(grid, u, u_src, ke_edge, h_edge, tend_u, err)
+ !$OMP SINGLE
call mpas_timer_stop("forcings", velForceTimer)
+ !$OMP END SINGLE
!
! velocity tendency: vertical mixing d/dz( nu_v du/dz))
!
+ !$OMP SINGLE
if (.not.config_implicit_vertical_mix) then
call mpas_timer_start("explicit vmix", .false., velExpVmixTimer)
call ocn_vel_vmix_tend_explicit(grid, u, h_edge, vertvisctopofedge, tend_u, err)
call mpas_timer_stop("explicit vmix", velExpVmixTimer)
endif
+ !$OMP END SINGLE
+ !$OMP SINGLE
call mpas_timer_stop("ocn_tend_u")
+ !$OMP END SINGLE
end subroutine ocn_tend_u!}}}
@@ -275,13 +309,14 @@
!> This routine computes the scalar (tracer) tendency for the ocean
!
!-----------------------------------------------------------------------
- subroutine ocn_tend_scalar(tend, s, d, grid, dt)!{{{
+ subroutine ocn_tend_scalar(tend, s, d, grid, scratch, dt)!{{{
implicit none
type (tend_type), intent(inout) :: tend !< Input/Output: Tendency structure
type (state_type), intent(in) :: s !< Input: State information
type (diagnostics_type), intent(in) :: d !< Input: Diagnostic information
type (mesh_type), intent(in) :: grid !< Input: Grid information
+ type (scratch_type), intent(inout) :: scratch !< Input/Output: Scratch structure
real (kind=RKIND), intent(in) :: dt !< Input: Time step
real (kind=RKIND), dimension(:,:), pointer :: &
@@ -291,7 +326,9 @@
integer :: err, iEdge, k
+ !$OMP SINGLE
call mpas_timer_start("ocn_tend_scalar")
+ !$OMP END SINGLE
uTransport => s % uTransport % array
h => s % h % array
@@ -303,20 +340,28 @@
tend_tr => tend % tracers % array
tend_h => tend % h % array
- allocate(uh(grid % nVertLevels, grid % nEdges+1))
+ !$OMP SINGLE
+ call mpas_allocate_scratch_field(scratch % uh)
+ !$OMP END SINGLE
+
+ uh => scratch % uh % array
!
! QC Comment (3/15/12): need to make sure that uTransport is the right
! transport velocity for the tracer.
+ !$OMP DO PRIVATE(k)
do iEdge = 1, grid % nEdges
do k = 1, grid % nVertLevels
uh(k, iEdge) = uTransport(k, iEdge) * h_edge(k, iEdge)
end do
end do
+ !$OMP END DO
!
! initialize tracer tendency (RHS of tracer equation) to zero.
!
+ !$OMP WORKSHARE
tend_tr(:,:,:) = 0.0
+ !$OMP END WORKSHARE
!
! tracer tendency: horizontal advection term -div( h \phi u)
@@ -327,16 +372,24 @@
! tracer_edge at the boundary will also need to be defined for flux boundaries.
! Monotonoic Advection, or standard advection
+ !$OMP SINGLE
call mpas_timer_start("adv", .false., tracerHadvTimer)
- call mpas_ocn_tracer_advection_tend(tracers, uh, wTop, h, h, dt, grid, tend_h, tend_tr)
+ !$OMP END SINGLE
+ call mpas_ocn_tracer_advection_tend(tracers, uh, wTop, h, h, dt, grid, tend_h, tend_tr, scratch)
+ !$OMP SINGLE
call mpas_timer_stop("adv", tracerHadvTimer)
+ !$OMP END SINGLE
!
! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="gray">abla \phi)
!
+ !$OMP SINGLE
call mpas_timer_start("hmix", .false., tracerHmixTimer)
- call ocn_tracer_hmix_tend(grid, h_edge, tracers, tend_tr, err)
+ !$OMP END SINGLE
+ call ocn_tracer_hmix_tend(grid, h_edge, tracers, tend_tr, scratch, err)
+ !$OMP SINGLE
call mpas_timer_stop("hmix", tracerHmixTimer)
+ !$OMP END SINGLE
! mrp 110516 printing
!print *, 'tend_tr 1',minval(tend_tr(3,1,1:nCells)),&
@@ -348,6 +401,7 @@
!
! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
!
+ !$OMP SINGLE
if (.not.config_implicit_vertical_mix) then
call mpas_timer_start("explicit vmix", .false., tracerExpVmixTimer)
@@ -355,6 +409,7 @@
call mpas_timer_stop("explicit vmix", tracerExpVmixTimer)
endif
+ !$OMP END SINGLE
! mrp 110516 printing
!print *, 'tend_tr 2',minval(tend_tr(3,1,1:nCells)),&
@@ -364,17 +419,24 @@
!
! add restoring to T and S in top model layer
!
+ !$OMP SINGLE
call mpas_timer_start("restoring", .false., tracerRestoringTimer)
+ !$OMP END SINGLE
call ocn_restoring_tend(grid, h, s%index_temperature, s%index_salinity, tracers, tend_tr, err)
+ !$OMP SINGLE
call mpas_timer_stop("restoring", tracerRestoringTimer)
- 10 format(2i8,10e20.10)
call mpas_timer_stop("ocn_tend_scalar")
+ !$OMP END SINGLE
- deallocate(uh)
+ !$OMP SINGLE
+ call mpas_deallocate_scratch_field(scratch % uh)
+ !$OMP END SINGLE
+ 10 format(2i8,10e20.10)
+
end subroutine ocn_tend_scalar!}}}
!***********************************************************************
@@ -491,9 +553,12 @@
! assign h_edge for maxLevelEdgeTop:maxLevelEdgeBot in the following section
! initialize h_edge to avoid divide by zero and NaN problems.
+ !$OMP WORKSHARE
h_edge = -1.0e34
+ !$OMP END WORKSHARE
coef_3rd_order = config_coef_3rd_order
+ !$OMP DO PRIVATE(cell1, cell2, k)
do iEdge=1,nEdges*hadv2nd
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -501,7 +566,9 @@
h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
end do
end do
+ !$OMP END DO
+ !$OMP DO PRIVATE(cell1, cell2, k, d2fdx2_cell1, d2fdx2_cell2, boundaryMask, i, velMask)
do iEdge=1,nEdges*hadv3rd
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -535,7 +602,9 @@
end do ! do k
end do ! do iEdge
+ !$OMP END DO
+ !$OMP DO PRIVATE(cell1, cell2, k, d2fdx2_cell1, d2fdx2_cell2, boundaryMask, i)
do iEdge=1,nEdges*hadv4th
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -568,11 +637,13 @@
end do ! do k
end do ! do iEdge
+ !$OMP END DO
!
! set the velocity and height at dummy address
! used -1e34 so error clearly occurs if these values are used.
!
+ !$OMP WORKSHARE
u(:,nEdges+1) = -1e34
h(:,nCells+1) = -1e34
tracers(s % index_temperature,:,nCells+1) = -1e34
@@ -583,6 +654,9 @@
divergence(:,:) = 0.0
ke(:,:) = 0.0
v(:,:) = 0.0
+ !$OMP END WORKSHARE
+
+ !$OMP DO PRIVATE(invAreaTri1, i, iEdge, k, r_tmp)
do iVertex = 1, nVertices
invAreaTri1 = 1.0 / areaTriangle(iVertex)
do i = 1, vertexDegree
@@ -595,7 +669,9 @@
end do
end do
end do
+ !$OMP END DO
+ !$OMP DO PRIVATE(invAreaCell1, i, iEdge, k, r_tmp)
do iCell = 1, nCells
invAreaCell1 = 1.0 / areaCell(iCell)
do i = 1, nEdgesOnCell(iCell)
@@ -608,7 +684,9 @@
end do
end do
end do
+ !$OMP END DO
+ !$OMP DO PRIVATE(i, eoe, k)
do iEdge=1,nEdges
! Compute v (tangential) velocities
do i=1,nEdgesOnEdge(iEdge)
@@ -619,13 +697,17 @@
v(k,iEdge) = v(k,iEdge) + weightsOnEdge(i,iEdge) * u(k, eoe)
end do
end do
-
end do
+ !$OMP END DO
!
! Compute kinetic energy in each vertex
!
+ !$OMP WORKSHARE
kev(:,:) = 0.0; kevc(:,:) = 0.0
+ !$OMP END WORKSHARE
+
+ !$OMP DO PRIVATE(i, iEdge, r_tmp, k)
do iVertex = 1, nVertices*ke_vertex_flag
do i = 1, vertexDegree
iEdge = edgesOnVertex(i, iVertex)
@@ -635,7 +717,9 @@
end do
end do
end do
+ !$OMP END DO
+ !$OMP DO PRIVATE(invAreaCell1, i, j, iVertex, k)
do iCell = 1, nCells*ke_vertex_flag
invAreaCell1 = 1.0 / areaCell(iCell)
do i = 1, nEdgesOnCell(iCell)
@@ -646,21 +730,25 @@
end do
end do
end do
+ !$OMP END DO
!
! Compute kinetic energy in each cell by blending ke and kevc
!
+ !$OMP DO PRIVATE(k)
do iCell=1,nCells*ke_vertex_flag
do k=1,nVertLevels
ke(k,iCell) = 5.0/8.0*ke(k,iCell) + 3.0/8.0*kevc(k,iCell)
end do
end do
+ !$OMP END DO
!
! Compute ke on cell edges at velocity locations for quadratic bottom drag.
!
! mrp 101025 efficiency note: we could get rid of ke_edge completely by
! using sqrt(u(k,iEdge)**2 + v(k,iEdge)**2) in its place elsewhere.
+ !$OMP DO PRIVATE(cell1, cell2, k)
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -668,11 +756,13 @@
ke_edge(k,iEdge) = 0.5 * (ke(k,cell1) + ke(k,cell2))
end do
end do
+ !$OMP END DO
!
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes Vor_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
+ !$OMP DO PRIVATE(invAreaTri1, k, h_vertex, i)
do iVertex = 1,nVertices
invAreaTri1 = 1.0 / areaTriangle(iVertex)
do k=1,maxLevelVertexBot(iVertex)
@@ -685,9 +775,14 @@
Vor_vertex(k,iVertex) = (fCoef*fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
end do
end do
+ !$OMP END DO
+ !$OMP WORKSHARE
Vor_cell(:,:) = 0.0
Vor_edge(:,:) = 0.0
+ !$OMP END WORKSHARE
+
+ !$OMP DO PRIVATE(vertex1, vertex2, k)
do iEdge = 1, nEdges
vertex1 = verticesOnEdge(1, iEdge)
vertex2 = verticesOnEdge(2, iEdge)
@@ -695,7 +790,9 @@
Vor_edge(k, iEdge) = 0.5 * (Vor_vertex(k, vertex1) + Vor_vertex(k, vertex2))
end do
end do
+ !$OMP END DO
+ !$OMP DO PRIVATE(invAreaCell1, i, j, iVertex, k)
do iCell = 1, nCells
invAreaCell1 = 1.0 / areaCell(iCell)
@@ -707,7 +804,9 @@
end do
end do
end do
+ !$OMP END DO
+ !$OMP DO PRIVATE(cell1, cell2, vertex1, vertex2, invLength, k)
do iEdge = 1,nEdges
cell1 = cellsOnEdge(1, iEdge)
cell2 = cellsOnEdge(2, iEdge)
@@ -727,12 +826,13 @@
do k = 1,maxLevelEdgeBot(iEdge)
gradVor_t(k,iEdge) = (Vor_vertex(k,vertex2) - Vor_vertex(k,vertex1)) * invLength
enddo
-
enddo
+ !$OMP END DO
!
! Modify PV edge with upstream bias.
!
+ !$OMP DO PRIVATE(k)
do iEdge = 1,nEdges
do k = 1,maxLevelEdgeBot(iEdge)
Vor_edge(k,iEdge) = Vor_edge(k,iEdge) &
@@ -740,6 +840,7 @@
+ v(k,iEdge) * gradVor_t(k,iEdge) )
enddo
enddo
+ !$OMP END DO
!
! equation of state
@@ -747,11 +848,14 @@
! For an isopycnal model, density should remain constant.
! For zlevel, calculate in-situ density
if (config_vert_grid_type.ne.'isopycnal') then
+ !DWJ 01/29/13 OMP Equation Of State Call....
+ !$OMP SINGLE
call mpas_timer_start("equation of state", .false., diagEOSTimer)
call ocn_equation_of_state_rho(s, grid, 0, 'relative', err)
! mrp 110324 In order to visualize rhoDisplaced, include the following
call ocn_equation_of_state_rho(s, grid, 1, 'relative', err)
call mpas_timer_stop("equation of state", diagEOSTimer)
+ !$OMP END SINGLE
endif
!
@@ -765,6 +869,7 @@
! Compute pressure at top of each layer, and then
! Montgomery Potential.
allocate(pTop(nVertLevels))
+ !$OMP DO PRIVATE(k)
do iCell=1,nCells
! assume atmospheric pressure at the surface is zero for now.
@@ -784,10 +889,12 @@
end do
end do
+ !$OMP END DO
deallocate(pTop)
else
+ !$OMP DO PRIVATE(k)
do iCell=1,nCells
! pressure for generalized coordinates
! assume atmospheric pressure at the surface is zero for now.
@@ -814,12 +921,14 @@
end do
end do
+ !$OMP END DO
endif
!
! Sea Surface Height
!
+ !$OMP DO
do iCell=1,nCells
! Start at the bottom where we know the depth, and go up.
! The bottom depth for this cell is bottomDepth(iCell).
@@ -827,8 +936,8 @@
! and z-coordinates are negative below the surface.
ssh(iCell) = - bottomDepth(iCell) + sum(h(1:maxLevelCell(iCell),iCell))
-
end do
+ !$OMP END DO
!
! Apply the GM closure as a bolus velocity
@@ -837,7 +946,9 @@
call ocn_gm_compute_uBolus(s,grid)
else
! mrp efficiency note: if uBolusGM is guaranteed to be zero, this can be removed.
+ !$OMP WORKSHARE
uBolusGM = 0.0
+ !$OMP END WORKSHARE
end if
end subroutine ocn_diagnostic_solve!}}}
@@ -928,7 +1039,9 @@
if (config_vert_grid_type.eq.'isopycnal') then
! set vertical velocity to zero in isopycnal case
+ !$OMP WORKSHARE
wTop=0.0_RKIND
+ !$OMP END WORKSHARE
return
end if
@@ -939,6 +1052,7 @@
! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
!
+ !$OMP DO PRIVATE(hSum, invAreaCell, i, iEdge, k, flux)
do iCell=1,nCells
div_hu(:) = 0.0_RKIND
div_hu_btr = 0.0_RKIND
@@ -973,6 +1087,7 @@
wTop(k,iCell) = wTop(k+1,iCell) - div_hu(k) - h_tend_col(k)
end do
end do
+ !$OMP END DO
deallocate(div_hu, h_tend_col)
@@ -1009,7 +1124,9 @@
integer, dimension(:), pointer :: maxLevelEdgeTop, nEdgesOnEdge
integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnEdge
+ !$OMP SINGLE
call mpas_timer_start("ocn_fuperp")
+ !$OMP END SINGLE
u => s % u % array
uBcl => s % uBcl % array
@@ -1027,6 +1144,7 @@
!
! Put f*uBcl^{perp} in u as a work variable
!
+ !$OMP DO PRIVATE(cell1, cell2, k, j, eoe)
do iEdge=1,nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -1040,8 +1158,11 @@
end do
end do
end do
+ !$OMP END DO
+ !$OMP SINGLE
call mpas_timer_stop("ocn_fuperp")
+ !$OMP END SINGLE
end subroutine ocn_fuperp!}}}
@@ -1068,13 +1189,16 @@
real (kind=RKIND), dimension(:,:), pointer :: h_edge, u
integer, dimension(:), pointer :: maxLevelEdgeTop
+ !$OMP SINGLE
call mpas_timer_start("ocn_filter_btr_mode_u")
+ !$OMP END SINGLE
u => s % u % array
h_edge => s % h_edge % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
nEdges = grid % nEdges
+ !$OMP DO PRIVATE(uhSum, hSum, k, vertSum)
do iEdge=1,nEdges
! hSum is initialized outside the loop because on land boundaries
@@ -1093,8 +1217,11 @@
u(k,iEdge) = u(k,iEdge) - vertSum
enddo
enddo ! iEdge
+ !$OMP END DO
+ !$OMP SINGLE
call mpas_timer_stop("ocn_filter_btr_mode_u")
+ !$OMP END SINGLE
end subroutine ocn_filter_btr_mode_u!}}}
@@ -1123,13 +1250,16 @@
integer, dimension(:), pointer :: maxLevelEdgeTop
+ !$OMP SINGLE
call mpas_timer_start("ocn_filter_btr_mode_tend_u")
+ !$OMP END SINGLE
tend_u => tend % u % array
h_edge => s % h_edge % array
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
nEdges = grid % nEdges
+ !$OMP DO PRIVATE(uhSum, hSum, k, vertSum)
do iEdge=1,nEdges
! hSum is initialized outside the loop because on land boundaries
@@ -1148,8 +1278,11 @@
tend_u(k,iEdge) = tend_u(k,iEdge) - vertSum
enddo
enddo ! iEdge
+ !$OMP END DO
+ !$OMP SINGLE
call mpas_timer_stop("ocn_filter_btr_mode_tend_u")
+ !$OMP END SINGLE
end subroutine ocn_filter_btr_mode_tend_u!}}}
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_thick_hadv.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_thick_hadv.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_thick_hadv.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -134,6 +134,7 @@
edgesOnCell => grid % edgesOnCell % array
edgeSignOnCell => grid % edgeSignOnCell % array
+ !$OMP DO PRIVATE(invAreaCell, i, iEdge, k, flux)
do iCell = 1, nCells
invAreaCell = 1.0 / areaCell(iCell)
do i = 1, nEdgesOnCell(iCell)
@@ -144,6 +145,7 @@
end do
end do
end do
+ !$OMP END DO
!--------------------------------------------------------------------
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_thick_vadv.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_thick_vadv.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_thick_vadv.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -116,11 +116,13 @@
nCells = grid % nCells
nVertLevels = grid % nVertLevels
+ !$OMP DO PRIVATE(k)
do iCell=1,nCells
do k=1,maxLevelCell(iCell)
tend(k,iCell) = tend(k,iCell) + wTop(k+1,iCell) - wTop(k,iCell)
end do
end do
+ !$OMP END DO
!--------------------------------------------------------------------
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_time_integration.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_time_integration.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_time_integration.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -96,6 +96,7 @@
call ocn_time_integrator_split(domain, dt)
endif
+ !$OMP SINGLE
block => domain % blocklist
do while (associated(block))
block % state % time_levs(2) % state % xtime % scalar = timeStamp
@@ -109,6 +110,7 @@
block => block % next
end do
+ !$OMP END SINGLE
end subroutine ocn_timestep!}}}
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_time_integration_rk4.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -93,7 +93,9 @@
real (kind=RKIND), dimension(:), allocatable:: A,C,uTemp
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
+ !$OMP SINGLE
call mpas_setup_provis_states(domain % blocklist)
+ !$OMP END SINGLE
!
! Initialize time_levs(2) with state at current time
@@ -103,16 +105,23 @@
!
block => domain % blocklist
do while (associated(block))
+ !$OMP WORKSHARE
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ !$OMP END WORKSHARE
+
+ !$OMP DO PRIVATE(k)
do iCell=1,block % mesh % nCells ! couple tracers to h
do k=1,block % mesh % maxLevelCell % array(iCell)
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = block % state % time_levs(1) % state % tracers % array(:,k,iCell) &
* block % state % time_levs(1) % state % h % array(k,iCell)
end do
end do
+ !$OMP END DO
+ !$OMP SINGLE
call mpas_copy_state(block % provis, block % state % time_levs(1) % state)
+ !$OMP END SINGLE
block => block % next
end do
@@ -127,14 +136,16 @@
rk_substep_weights(3) = dt
rk_substep_weights(4) = 0.
-
+ !$OMP SINGLE
call mpas_timer_start("RK4-main loop")
+ !$OMP END SINGLE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN RK loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do rk_step = 1, 4
! --- update halos for diagnostic variables
+ !$OMP SINGLE
call mpas_timer_start("RK4-diagnostic halo update")
call mpas_dmpar_exch_halo_field(domain % blocklist % provis % Vor_edge)
if (config_h_mom_eddy_visc4 > 0.0) then
@@ -142,33 +153,37 @@
call mpas_dmpar_exch_halo_field(domain % blocklist % provis % vorticity)
end if
call mpas_timer_stop("RK4-diagnostic halo update")
+ !$OMP END SINGLE
! --- compute tendencies
+ !$OMP SINGLE
call mpas_timer_start("RK4-tendency computations")
+ !$OMP END SINGLE
block => domain % blocklist
do while (associated(block))
if (.not.config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, block % provis, block % diagnostics, err)
+ call ocn_vmix_coefs(block % mesh, block % provis, block % diagnostics, block % scratch, err)
end if
! advection of u uses u, while advection of h and tracers use uTransport.
call ocn_wtop(block % mesh, block % provis % h % array, block % provis % h_edge % array, &
block % provis % u % array, block % provis % wTop % array, err)
- call ocn_tend_u(block % tend, block % provis, block % diagnostics, block % mesh)
+ call ocn_tend_u(block % tend, block % provis, block % diagnostics, block % mesh, block % scratch)
call ocn_wtop(block % mesh, block % provis % h % array, block % provis % h_edge % array, &
block % provis % uTransport % array, block % provis % wTop % array, err)
- call ocn_tend_h(block % tend, block % provis, block % mesh)
+ call ocn_tend_h(block % tend, block % provis, block % mesh, block % scratch)
if (config_rk_filter_btr_mode) then
call ocn_filter_btr_mode_tend_u(block % tend, block % provis, block % mesh)
endif
- call ocn_tend_scalar(block % tend, block % provis, block % diagnostics, block % mesh, dt)
+ call ocn_tend_scalar(block % tend, block % provis, block % diagnostics, block % mesh, block % scratch, dt)
block => block % next
end do
+ !$OMP SINGLE
call mpas_timer_stop("RK4-tendency computations")
! --- update halos for prognostic variables
@@ -182,15 +197,20 @@
! --- compute next substep state
call mpas_timer_start("RK4-update diagnostic variables")
+ !$OMP END SINGLE
if (rk_step < 4) then
block => domain % blocklist
do while (associated(block))
+ !$OMP WORKSHARE
block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:) &
+ rk_substep_weights(rk_step) * block % tend % u % array(:,:)
block % provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:) &
+ rk_substep_weights(rk_step) * block % tend % h % array(:,:)
+ !$OMP END WORKSHARE
+
+ !$OMP DO PRIVATE(k)
do iCell=1,block % mesh % nCells
do k=1,block % mesh % maxLevelCell % array(iCell)
block % provis % tracers % array(:,k,iCell) = ( block % state % time_levs(1) % state % h % array(k,iCell) * &
@@ -200,41 +220,58 @@
end do
end do
+ !$OMP END DO
+
if (config_test_case == 1) then ! For case 1, wind field should be fixed
+ !$OMP WORKSHARE
block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ !$OMP END WORKSHARE
end if
if (config_prescribe_velocity) then
+ !$OMP WORKSHARE
block % provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ !$OMP END WORKSHARE
end if
if (config_prescribe_thickness) then
+ !$OMP WORKSHARE
block % provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ !$OMP END WORKSHARE
end if
call ocn_diagnostic_solve(dt, block % provis, block % mesh)
! Compute velocity transport, used in advection terms of h and tracer tendency
+ !$OMP WORKSHARE
block % provis % uTransport % array(:,:) &
= block % provis % u % array(:,:) &
+ block % provis % uBolusGM % array(:,:)
+ !$OMP END WORKSHARE
block => block % next
end do
end if
+ !$OMP SINGLE
call mpas_timer_stop("RK4-update diagnostic variables")
+ !$OMP END SINGLE
!--- accumulate update (for RK4)
+ !$OMP SINGLE
call mpas_timer_start("RK4-RK4 accumulate update")
+ !$OMP END SINGLE
block => domain % blocklist
do while (associated(block))
+ !$OMP WORKSHARE
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(2) % state % u % array(:,:) &
+ rk_weights(rk_step) * block % tend % u % array(:,:)
block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(2) % state % h % array(:,:) &
+ rk_weights(rk_step) * block % tend % h % array(:,:)
+ !$OMP END WORKSHARE
+ !$OMP DO PRIVATE(k)
do iCell=1,block % mesh % nCells
do k=1,block % mesh % maxLevelCell % array(iCell)
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
@@ -242,37 +279,48 @@
+ rk_weights(rk_step) * block % tend % tracers % array(:,k,iCell)
end do
end do
+ !$OMP END DO
block => block % next
end do
+ !$OMP SINGLE
call mpas_timer_stop("RK4-RK4 accumulate update")
+ !$OMP END SINGLE
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END RK loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !$OMP SINGLE
call mpas_timer_stop("RK4-main loop")
+ !$OMP END SINGLE
!
! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
!
+ !$OMP SINGLE
call mpas_timer_start("RK4-cleaup phase")
+ !$OMP END SINGLE
! Rescale tracers
block => domain % blocklist
do while(associated(block))
+ !$OMP DO PRIVATE(k)
do iCell = 1, block % mesh % nCells
do k = 1, block % mesh % maxLevelCell % array(iCell)
block % state % time_levs(2) % state % tracers % array(:, k, iCell) = block % state % time_levs(2) % state % tracers % array(:, k, iCell) &
/ block % state % time_levs(2) % state % h % array(k, iCell)
end do
end do
+ !$OMP END DO
block => block % next
end do
if (config_implicit_vertical_mix) then
+ !$OMP SINGLE
call mpas_timer_start("RK4-implicit vert mix")
+ !$OMP END SINGLE
block => domain % blocklist
do while(associated(block))
@@ -284,7 +332,7 @@
! implicit vmix routine that follows. mrp 121023.
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
- call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+ call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, block % scratch, err)
block => block % next
end do
@@ -292,35 +340,47 @@
! this leads to lack of volume conservation. It is required because halo updates in RK4 are only
! conducted on tendencies, not on the velocity and tracer fields. So this update is required to
! communicate the change due to implicit vertical mixing across the boundary.
+ !$OMP SINGLE
call mpas_timer_start("RK4-implicit vert mix halos")
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % u)
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % tracers)
call mpas_timer_stop("RK4-implicit vert mix halos")
call mpas_timer_stop("RK4-implicit vert mix")
+ !$OMP END SINGLE
end if
block => domain % blocklist
do while (associated(block))
if (config_test_case == 1) then ! For case 1, wind field should be fixed
+ !$OMP WORKSHARE
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ !$OMP END WORKSHARE
end if
if (config_prescribe_velocity) then
+ !$OMP WORKSHARE
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ !$OMP END WORKSHARE
end if
if (config_prescribe_thickness) then
+ !$OMP WORKSHARE
block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ !$OMP END WORKSHARE
end if
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
! Compute velocity transport, used in advection terms of h and tracer tendency
+ !$OMP WORKSHARE
block % state % time_levs(2) % state % uTransport % array(:,:) &
= block % state % time_levs(2) % state % u % array(:,:) &
+ block % state % time_levs(2) % state % uBolusGM % array(:,:)
+ !$OMP END WORKSHARE
+ !$OMP SECTIONS
+ !$OMP SECTION
call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
block % state % time_levs(2) % state % uReconstructX % array, &
block % state % time_levs(2) % state % uReconstructY % array, &
@@ -330,6 +390,7 @@
)
!TDR
+ !$OMP SECTION
call mpas_reconstruct(block % mesh, block % mesh % u_src % array, &
block % state % time_levs(2) % state % uSrcReconstructX % array, &
block % state % time_levs(2) % state % uSrcReconstructY % array, &
@@ -337,15 +398,18 @@
block % state % time_levs(2) % state % uSrcReconstructZonal % array, &
block % state % time_levs(2) % state % uSrcReconstructMeridional % array &
)
+ !$OMP END SECTIONS
!TDR
call ocn_time_average_accumulate(block % state % time_levs(2) % state, block % state % time_levs(1) % state)
block => block % next
end do
+ !$OMP SINGLE
call mpas_timer_stop("RK4-cleaup phase")
call mpas_deallocate_provis_states(domain % blocklist)
+ !$OMP END SINGLE
end subroutine ocn_time_integrator_rk4!}}}
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_time_integration_split.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_time_integration_split.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -99,6 +99,7 @@
real (kind=RKIND), dimension(:), allocatable:: uTemp
real (kind=RKIND), dimension(:,:), allocatable:: tracersTemp
+ !$OMP SINGLE
call mpas_timer_start("se timestep", .false., timer_main)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -107,10 +108,12 @@
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call mpas_timer_start("se prep", .false., timer_prep)
+ !$OMP END SINGLE
block => domain % blocklist
do while (associated(block))
! Initialize * variables that are used to compute baroclinic tendencies below.
+ !$OMP DO PRIVATE(k)
do iEdge=1,block % mesh % nEdges
do k=1,block % mesh % nVertLevels !maxLevelEdgeTop % array(iEdge)
@@ -135,10 +138,14 @@
end do
end do
+ !$OMP END DO
+ !$OMP WORKSHARE
block % state % time_levs(2) % state % ssh % array(:) &
= block % state % time_levs(1) % state % ssh % array(:)
+ !$OMP END WORKSHARE
+ !$OMP DO PRIVATE(k)
do iCell=1,block % mesh % nCells
do k=1,block % mesh % maxLevelCell % array(iCell)
@@ -150,11 +157,14 @@
end do
end do
+ !$OMP END DO
block => block % next
end do
+ !$OMP SINGLE
call mpas_timer_stop("se prep", timer_prep)
+ !$OMP END SINGLE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN large iteration loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -165,6 +175,7 @@
do split_explicit_step = 1, config_n_ts_iter
! --- update halos for diagnostic variables
+ !$OMP SINGLE
call mpas_timer_start("se halo diag", .false., timer_halo_diagnostic)
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % Vor_edge)
if (config_h_mom_eddy_visc4 > 0.0) then
@@ -172,6 +183,7 @@
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % vorticity)
end if
call mpas_timer_stop("se halo diag", timer_halo_diagnostic)
+ !$OMP END SINGLE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
@@ -180,16 +192,20 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! compute velocity tendencies, T(u*,w*,p*)
+ !$OMP SINGLE
call mpas_timer_start("se bcl vel", .false., timer_bcl_vel)
+ !$OMP END SINGLE
block => domain % blocklist
do while (associated(block))
stage1_tend_time = min(split_explicit_step,2)
+ !$OMP SINGLE
if (.not.config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, block % state % time_levs(stage1_tend_time) % state, block % diagnostics, err)
+ call ocn_vmix_coefs(block % mesh, block % state % time_levs(stage1_tend_time) % state, block % diagnostics, block % scratch, err)
end if
+ !$OMP END SINGLE
! compute wTop. Use u (rather than uTransport) for momentum advection.
! Use the most recent time level available.
@@ -198,7 +214,7 @@
block % state % time_levs(stage1_tend_time) % state % u % array, &
block % state % time_levs(stage1_tend_time) % state % wTop % array, err)
- call ocn_tend_u(block % tend, block % state % time_levs(stage1_tend_time) % state, block % diagnostics, block % mesh)
+ call ocn_tend_u(block % tend, block % state % time_levs(stage1_tend_time) % state, block % diagnostics, block % mesh, block % scratch)
block => block % next
end do
@@ -222,6 +238,7 @@
! Put f*uBcl^{perp} in uNew as a work variable
call ocn_fuperp(block % state % time_levs(2) % state , block % mesh)
+ !$OMP DO PRIVATE(cell1, cell2, k, uhSum, hSum)
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -265,18 +282,22 @@
enddo
enddo ! iEdge
+ !$OMP END DO
deallocate(uTemp)
block => block % next
end do
+ !$OMP SINGLE
call mpas_timer_start("se halo ubcl", .false., timer_halo_ubcl)
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % uBcl)
call mpas_timer_stop("se halo ubcl", timer_halo_ubcl)
+ !$OMP END SINGLE
end do ! do j=1,config_n_bcl_iter
+ !$OMP SINGLE
call mpas_timer_stop("se bcl vel", timer_bcl_vel)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END baroclinic iterations on linear Coriolis term
@@ -290,6 +311,7 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
call mpas_timer_start("se btr vel", .false., timer_btr_vel)
+ !$OMP END SINGLE
oldBtrSubcycleTime = 1
newBtrSubcycleTime = 2
@@ -300,10 +322,13 @@
do while (associated(block))
! For Split_Explicit unsplit, simply set uBtrNew=0, uBtrSubcycle=0, and uNew=uBclNew
+ !$OMP WORKSHARE
block % state % time_levs(2) % state % uBtr % array(:) = 0.0
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(2) % state % uBcl % array(:,:)
+ !$OMP END WORKSHARE
+ !$OMP DO PRIVATE(k)
do iEdge=1,block % mesh % nEdges
do k=1,block % mesh % nVertLevels
@@ -317,6 +342,7 @@
enddo
end do ! iEdge
+ !$OMP END DO
block => block % next
end do ! block
@@ -328,15 +354,20 @@
do while (associated(block))
if (config_filter_btr_mode) then
+ !$OMP WORKSHARE
block % state % time_levs(1) % state % GBtrForcing % array(:) = 0.0
+ !$OMP END WORKSHARE
endif
+ !$OMP DO
do iCell=1,block % mesh % nCells
! sshSubcycleOld = sshOld
block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
= block % state % time_levs(1) % state % ssh % array(iCell)
end do
+ !$OMP END DO
+ !$OMP DO
do iEdge=1,block % mesh % nEdges
! uBtrSubcycleOld = uBtrOld
@@ -350,6 +381,7 @@
! FBtr = 0
block % state % time_levs(1) % state % FBtr % array(iEdge) = 0.0
end do
+ !$OMP END DO
block => block % next
end do ! block
@@ -362,12 +394,14 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Barotropic subcycle: VELOCITY PREDICTOR STEP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !$OMP BARRIER
if (config_btr_gam1_uWt1>1.0e-12) then ! only do this part if it is needed in next SSH solve
uPerpTime = oldBtrSubcycleTime
block => domain % blocklist
do while (associated(block))
+ !$OMP DO PRIVATE(cell1, cell2, CoriolisTerm, i, eoe)
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
@@ -392,14 +426,17 @@
/ block % mesh % dcEdge % array(iEdge) &
+ block % state % time_levs(1) % state % GBtrForcing % array(iEdge))) * block % mesh % edgeMask % array(1, iEdge)
end do
+ !$OMP END DO
block => block % next
end do ! block
! boundary update on uBtrNew
+ !$OMP SINGLE
call mpas_timer_start("se halo ubtr", .false., timer_halo_ubtr)
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle)
call mpas_timer_stop("se halo ubtr", timer_halo_ubtr)
+ !$OMP END SINGLE
endif ! config_btr_gam1_uWt1>1.0e-12
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -408,7 +445,9 @@
block => domain % blocklist
do while (associated(block))
+ !$OMP WORKSHARE
block % tend % ssh % array(:) = 0.0
+ !$OMP END WORKSHARE
if (config_btr_solve_SSH2) then
! If config_btr_solve_SSH2=.true., then do NOT accumulate FBtr in this SSH predictor
@@ -425,6 +464,7 @@
! config_btr_gam1_uWt1= 0 flux = uBtrOld*H
! mrp 120201 efficiency: could we combine the following edge and cell loops?
+ !$OMP DO PRIVATE(i, iEdge, cell1, cell2, sshEdge, hSum, flux)
do iCell = 1, block % mesh % nCells
do i = 1, block % mesh % nEdgesOnCell % array(iCell)
iEdge = block % mesh % edgesOnCell % array(i, iCell)
@@ -457,7 +497,9 @@
end do
end do
+ !$OMP END DO
+ !$OMP DO PRIVATE(cell1, cell2, sshEdge, hSum, flux)
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -484,8 +526,10 @@
block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) &
+ FBtr_coeff*flux
end do
+ !$OMP END DO
! SSHnew = SSHold + dt/J*(-div(Flux))
+ !$OMP DO PRIVATE(iCell)
do iCell=1,block % mesh % nCells
block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
@@ -493,25 +537,35 @@
+ dt/config_n_btr_subcycles * block % tend % ssh % array(iCell) / block % mesh % areaCell % array (iCell)
end do
+ !$OMP END DO
block => block % next
end do ! block
! boundary update on SSHnew
+ !$OMP SINGLE
call mpas_timer_start("se halo ssh", .false., timer_halo_ssh)
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle)
call mpas_timer_stop("se halo ssh", timer_halo_ssh)
+ !$OMP END SINGLE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Barotropic subcycle: VELOCITY CORRECTOR STEP
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do BtrCorIter=1,config_n_btr_cor_iter
uPerpTime = newBtrSubcycleTime
+
+ !$OMP SINGLE
+ call mpas_allocate_scratch_field(domain % blocklist % scratch % uTempEdges)
+ !$OMP END SINGLE
block => domain % blocklist
do while (associated(block))
- allocate(utemp(block % mesh % nEdges+1))
- uTemp(:) = block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:)
+ !$OMP WORKSHARE
+ block % scratch % uTempEdges % array(:) = block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(:)
+ !$OMP END WORKSHARE
+
+ !$OMP DO PRIVATE(cell1, cell2, CoriolisTerm, i, eoe, sshCell1, sshCell2)
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -521,8 +575,9 @@
do i = 1,block % mesh % nEdgesOnEdge % array(iEdge)
eoe = block % mesh % edgesOnEdge % array(i,iEdge)
CoriolisTerm = CoriolisTerm + block % mesh % weightsOnEdge % array(i,iEdge) &
+ !DWJ 01/29/13: OpenMP fixes
!* block % state % time_levs(uPerpTime) % state % uBtrSubcycle % array(eoe) &
- * uTemp(eoe) &
+ * block % scratch % uTempEdges % array(eoe) &
* block % mesh % fEdge % array(eoe)
end do
@@ -539,15 +594,21 @@
+ dt/config_n_btr_subcycles *(CoriolisTerm - gravity *(sshCell2 - sshCell1) /block % mesh % dcEdge % array(iEdge) &
+ block % state % time_levs(1) % state % GBtrForcing % array(iEdge))) * block % mesh % edgeMask % array(1,iEdge)
end do
- deallocate(uTemp)
-
+ !$OMP END DO
+
block => block % next
end do ! block
+
+ !$OMP SINGLE
+ call mpas_deallocate_scratch_field(domain % blocklist % scratch % uTempEdges)
+ !$OMP END SINGLE
! boundary update on uBtrNew
+ !$OMP SINGLE
call mpas_timer_start("se halo ubtr", .false., timer_halo_ubtr)
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle)
call mpas_timer_stop("se halo ubtr", timer_halo_ubtr)
+ !$OMP END SINGLE
end do !do BtrCorIter=1,config_n_btr_cor_iter
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -557,7 +618,9 @@
block => domain % blocklist
do while (associated(block))
+ !$OMP WORKSHARE
block % tend % ssh % array(:) = 0.0
+ !$OMP END WORKSHARE
! config_btr_gam3_uWt2 sets the forward weighting of velocity in the SSH computation
! config_btr_gam3_uWt2= 1 flux = uBtrNew*H
@@ -565,6 +628,7 @@
! config_btr_gam3_uWt2= 0 flux = uBtrOld*H
! mrp 120201 efficiency: could we combine the following edge and cell loops?
+ !$OMP DO PRIVATE(i, iEdge, cell1, cell2, sshCell1, sshCell2, sshEdge, hSum, flux)
do iCell = 1, block % mesh % nCells
do i = 1, block % mesh % nEdgesOnCell % array(iCell)
iEdge = block % mesh % edgesOnCell % array(i, iCell)
@@ -602,7 +666,9 @@
end do
end do
+ !$OMP END DO
+ !$OMP DO PRIVATE(cell1, cell2, sshCell1, sshCell2, sshEdge, hSum, flux)
do iEdge=1,block % mesh % nEdges
cell1 = block % mesh % cellsOnEdge % array(1,iEdge)
cell2 = block % mesh % cellsOnEdge % array(2,iEdge)
@@ -632,21 +698,26 @@
block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) + flux
end do
+ !$OMP END DO
! SSHnew = SSHold + dt/J*(-div(Flux))
+ !$OMP DO
do iCell=1,block % mesh % nCells
block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
= block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(iCell) &
+ dt/config_n_btr_subcycles * block % tend % ssh % array(iCell) / block % mesh % areaCell % array (iCell)
end do
+ !$OMP END DO
block => block % next
end do ! block
! boundary update on SSHnew
+ !$OMP SINGLE
call mpas_timer_start("se halo ssh", .false., timer_halo_ssh)
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle)
call mpas_timer_stop("se halo ssh", timer_halo_ssh)
+ !$OMP END SINGLE
endif ! config_btr_solve_SSH2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -660,6 +731,7 @@
! This accumulates the sum.
! If the Barotropic Coriolis iteration is limited to one, this could
! be merged with the above code.
+ !$OMP DO
do iEdge=1,block % mesh % nEdges
block % state % time_levs(2) % state % uBtr % array(iEdge) &
@@ -667,12 +739,14 @@
+ block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)
end do ! iEdge
+ !$OMP END DO
block => block % next
end do ! block
! advance time pointers
oldBtrSubcycleTime = mod(oldBtrSubcycleTime,2)+1
newBtrSubcycleTime = mod(newBtrSubcycleTime,2)+1
+ !$OMP BARRIER
end do ! j=1,config_n_btr_subcycles
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -683,6 +757,7 @@
block => domain % blocklist
do while (associated(block))
+ !$OMP DO
do iEdge=1,block % mesh % nEdges
block % state % time_levs(1) % state % FBtr % array(iEdge) = block % state % time_levs(1) % state % FBtr % array(iEdge) &
/ (config_n_btr_subcycles*config_btr_subcycle_loop_factor)
@@ -690,17 +765,20 @@
block % state % time_levs(2) % state % uBtr % array(iEdge) = block % state % time_levs(2) % state % uBtr % array(iEdge) &
/ (config_n_btr_subcycles*config_btr_subcycle_loop_factor + 1)
end do
+ !$OMP END DO
block => block % next
end do ! block
! boundary update on F
+ !$OMP BARRIER
+ !$OMP SINGLE
call mpas_timer_start("se halo F", .false., timer_halo_f)
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(1) % state % FBtr)
call mpas_timer_stop("se halo F", timer_halo_f)
+ !$OMP END SINGLE
-
! Check that you can compute SSH using the total sum or the individual increments
! over the barotropic subcycles.
! efficiency: This next block of code is really a check for debugging, and can
@@ -722,6 +800,7 @@
ucorr_coef = 0
endif
+ !$OMP DO PRIVATE(uhSum, hSum, k, uCorr)
do iEdge=1,block % mesh % nEdges
! velocity for uCorrection is uBtr + uBcl + uBolus
@@ -758,6 +837,7 @@
enddo
end do ! iEdge
+ !$OMP END DO
deallocate(uTemp)
@@ -766,7 +846,9 @@
endif ! split_explicit
+ !$OMP SINGLE
call mpas_timer_stop("se btr vel", timer_btr_vel)
+ !$OMP END SINGLE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
@@ -791,26 +873,30 @@
block % state % time_levs(2) % state % uTransport % array, &
block % state % time_levs(2) % state % wTop % array, err)
- call ocn_tend_h(block % tend, block % state % time_levs(2) % state, block % mesh)
+ call ocn_tend_h(block % tend, block % state % time_levs(2) % state, block % mesh, block % scratch)
block => block % next
end do
! update halo for thickness and tracer tendencies
+ !$OMP SINGLE
call mpas_timer_start("se halo h", .false., timer_halo_h)
call mpas_dmpar_exch_halo_field(domain % blocklist % tend % h)
call mpas_timer_stop("se halo h", timer_halo_h)
+ !$OMP END SINGLE
block => domain % blocklist
do while (associated(block))
- call ocn_tend_scalar(block % tend, block % state % time_levs(2) % state, block % diagnostics, block % mesh, dt)
+ call ocn_tend_scalar(block % tend, block % state % time_levs(2) % state, block % diagnostics, block % mesh, block % scratch, dt)
block => block % next
end do
! update halo for thickness and tracer tendencies
+ !$OMP SINGLE
call mpas_timer_start("se halo tracers", .false., timer_halo_tracers)
call mpas_dmpar_exch_halo_field(domain % blocklist % tend % tracers)
call mpas_timer_stop("se halo tracers", timer_halo_tracers)
+ !$OMP END SINGLE
block => domain % blocklist
do while (associated(block))
@@ -827,6 +913,7 @@
! Only need T & S for earlier iterations,
! then all the tracers needed the last time through.
+ !$OMP DO PRIVATE(k, temp_h, i)
do iCell=1,block % mesh % nCells
! sshNew is a pointer, defined above.
do k=1,block % mesh % maxLevelCell % array(iCell)
@@ -858,9 +945,10 @@
end do
end do
end do ! iCell
+ !$OMP END DO
+ !$OMP DO PRIVATE(k)
do iEdge=1,block % mesh % nEdges
-
do k=1,block % mesh % nVertLevels
! u = uBtr + uBcl
@@ -872,8 +960,8 @@
+ block % state % time_levs(2) % state % uBcl % array(k,iEdge) )
enddo
-
end do ! iEdge
+ !$OMP END DO
! mrp 110512 I really only need this to compute h_edge, density, pressure, and SSH
! I can par this down later.
@@ -889,6 +977,7 @@
!TDR: should we move this code into a subroutine called "compute_final_values_at_nplus1"?
!TDR: this could be within a contains statement in this routine
+ !$OMP DO PRIVATE(k, i)
do iCell=1,block % mesh % nCells
do k=1,block % mesh % maxLevelCell % array(iCell)
@@ -908,6 +997,7 @@
enddo
end do
end do
+ !$OMP END DO
! Recompute final u to go on to next step.
! u_{n+1} = uBtr_{n+1} + uBcl_{n+1}
@@ -918,6 +1008,7 @@
! note that uBcl is recomputed at the beginning of the next timestep due to Imp Vert mixing,
! so uBcl does not have to be recomputed here.
+ !$OMP DO PRIVATE(k)
do iEdge=1,block % mesh % nEdges
do k=1,block % mesh % maxLevelEdgeTop % array(iEdge)
block % state % time_levs(2) % state % u % array(k,iEdge) &
@@ -926,21 +1017,21 @@
- block % state % time_levs(1) % state % uBcl % array(k,iEdge)
end do
end do ! iEdges
+ !$OMP END DO
endif ! split_explicit_step
block => block % next
end do
-
-
-
end do ! split_explicit_step = 1, config_n_ts_iter
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END large iteration loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (config_implicit_vertical_mix) then
+ !$OMP SINGLE
call mpas_timer_start("se implicit vert mix")
+ !$OMP END SINGLE
block => domain % blocklist
do while(associated(block))
@@ -952,7 +1043,7 @@
! implicit vmix routine that follows. mrp 121023.
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
- call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+ call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, block % scratch, err)
block => block % next
end do
@@ -961,36 +1052,47 @@
! this leads to lack of volume conservation. It is required because halo updates in stage 3 are only
! conducted on tendencies, not on the velocity and tracer fields. So this update is required to
! communicate the change due to implicit vertical mixing across the boundary.
+ !$OMP SINGLE
call mpas_timer_start("se implicit vert mix halos")
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % u)
call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % tracers)
call mpas_timer_stop("se implicit vert mix halos")
call mpas_timer_stop("se implicit vert mix")
+ !$OMP END SINGLE
end if
block => domain % blocklist
do while (associated(block))
if (config_test_case == 1) then ! For case 1, wind field should be fixed
+ !$OMP WORKSHARE
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ !$OMP END WORKSHARE
end if
if (config_prescribe_velocity) then
+ !$OMP WORKSHARE
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ !$OMP END WORKSHARE
end if
if (config_prescribe_thickness) then
+ !$OMP WORKSHARE
block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ !$OMP END WORKSHARE
end if
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
! Compute velocity transport, used in advection terms of h and tracer tendency
+ !$OMP WORKSHARE
block % state % time_levs(2) % state % uTransport % array(:,:) &
= block % state % time_levs(2) % state % u % array(:,:) &
+ block % state % time_levs(2) % state % uBolusGM % array(:,:)
+ !$OMP END WORKSHARE
+ !$OMP SINGLE
call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
block % state % time_levs(2) % state % uReconstructX % array, &
block % state % time_levs(2) % state % uReconstructY % array, &
@@ -998,8 +1100,10 @@
block % state % time_levs(2) % state % uReconstructZonal % array, &
block % state % time_levs(2) % state % uReconstructMeridional % array &
)
+ !$OMP END SINGLE NOWAIT
!TDR
+ !$OMP SINGLE
call mpas_reconstruct(block % mesh, block % mesh % u_src % array, &
block % state % time_levs(2) % state % uSrcReconstructX % array, &
block % state % time_levs(2) % state % uSrcReconstructY % array, &
@@ -1007,6 +1111,7 @@
block % state % time_levs(2) % state % uSrcReconstructZonal % array, &
block % state % time_levs(2) % state % uSrcReconstructMeridional % array &
)
+ !$OMP END SINGLE
!TDR
call ocn_time_average_accumulate(block % state % time_levs(2) % state, block % state % time_levs(1) % state)
@@ -1014,7 +1119,9 @@
block => block % next
end do
+ !$OMP SINGLE
call mpas_timer_stop("se timestep", timer_main)
+ !$OMP END SINGLE
end subroutine ocn_time_integrator_split!}}}
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -239,7 +239,7 @@
!> advection of tracers.
!
!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend)!{{{
+ subroutine mpas_ocn_tracer_advection_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend, scratch)!{{{
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: tracer tendency
real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !< Input/Output: tracer values
@@ -250,9 +250,10 @@
real (kind=RKIND), intent(in) :: dt !< Input: Time step
type (mesh_type), intent(in) :: grid !< Input: grid information
real (kind=RKIND), dimension(:,:), intent(in) :: tend_h !< Input: Thickness tendency information
+ type (scratch_type), intent(inout) :: scratch !< Input/Output : Scratch structure
if(monotonicOn) then
- call mpas_ocn_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend)
+ call mpas_ocn_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend, scratch)
else
call mpas_ocn_tracer_advection_std_tend(tracers, uh, w, verticalCellSize, grid, tend)
endif
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -43,7 +43,7 @@
!> Both horizontal and vertical.
!
!-----------------------------------------------------------------------
- subroutine mpas_ocn_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend)!{{{
+ subroutine mpas_ocn_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend, scratch)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -60,6 +60,7 @@
real (kind=RKIND), intent(in) :: dt !< Input: Timestep
type (mesh_type), intent(in) :: grid !< Input: Grid information
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: Tracer tendency
+ type (scratch_type), intent(inout) :: scratch !< Input/Output: Scratch structure
integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
integer :: nVertLevels, num_tracers, nCells, nEdges, nCellsSolve
@@ -77,9 +78,41 @@
real (kind=RKIND), parameter :: eps = 1.e-10
- type (field2dReal), pointer :: high_order_horiz_flux_field
+ ! Initialize pointers
+ !$OMP SECTIONS
+ !$OMP SECTION
+ call mpas_allocate_scratch_field(scratch % high_order_horiz_flux)
+ !$OMP SECTION
+ call mpas_allocate_scratch_field(scratch % high_order_vert_flux)
+ !$OMP SECTION
+ call mpas_allocate_scratch_field(scratch % tracer_new)
+ !$OMP SECTION
+ call mpas_allocate_scratch_field(scratch % tracer_cur)
+ !$OMP SECTION
+ call mpas_allocate_scratch_field(scratch % upwind_tendency)
+ !$OMP SECTION
+ call mpas_allocate_scratch_field(scratch % inv_h_new)
+ !$OMP SECTION
+ call mpas_allocate_scratch_field(scratch % tracer_max)
+ !$OMP SECTION
+ call mpas_allocate_scratch_field(scratch % tracer_min)
+ !$OMP SECTION
+ call mpas_allocate_scratch_field(scratch % flux_incoming)
+ !$OMP SECTION
+ call mpas_allocate_scratch_field(scratch % flux_outgoing)
+ !$OMP END SECTIONS
- ! Initialize pointers
+ high_order_horiz_flux => scratch % high_order_horiz_flux % array
+ high_order_vert_flux => scratch % high_order_vert_flux % array
+ tracer_new => scratch % tracer_new % array
+ tracer_cur => scratch % tracer_cur % array
+ upwind_tendency => scratch % upwind_tendency % array
+ inv_h_new => scratch % inv_h_new % array
+ tracer_max => scratch % tracer_max % array
+ tracer_min => scratch % tracer_min % array
+ flux_incoming => scratch % flux_incoming % array
+ flux_outgoing => scratch % flux_outgoing % array
+
dvEdge => grid % dvEdge % array
cellsOnEdge => grid % cellsOnEdge % array
cellsOnCell => grid % cellsOnCell % array
@@ -104,43 +137,19 @@
nVertLevels = grid % nVertLevels
num_tracers = size(tracers,dim=1)
- allocate(high_order_horiz_flux_field)
- nullify(high_order_horiz_flux_field % next)
- high_order_horiz_flux_field % block => grid % block
- high_order_horiz_flux_field % sendList => grid % xEdge % sendList
- high_order_horiz_flux_field % recvList => grid % xEdge % recvList
- high_order_horiz_flux_field % copyList => grid % xEdge % copyList
- high_order_horiz_flux_field % dimSizes(1) = nVertLevels
- high_order_horiz_flux_field % dimSizes(2) = nEdges+1
- allocate(high_order_horiz_flux_field % array(high_order_horiz_flux_field % dimSizes(1), high_order_horiz_flux_field % dimSizes(2)))
- high_order_horiz_flux => high_order_horiz_flux_field % array
-
- ! allocate nCells arrays
-
- allocate(tracer_new(nVertLevels, nCells+1))
- allocate(tracer_cur(nVertLevels, nCells+1))
- allocate(upwind_tendency(nVertLevels, nCells+1))
- allocate(inv_h_new(nVertLevels, nCells+1))
- allocate(tracer_max(nVertLevels, nCells+1))
- allocate(tracer_min(nVertLevels, nCells+1))
- allocate(flux_incoming(nVertLevels, nCells+1))
- allocate(flux_outgoing(nVertLevels, nCells+1))
-
- ! allocate nEdges arrays
-! allocate(high_order_horiz_flux(nVertLevels, nEdges))
-
- ! allocate nVertLevels+1 and nCells arrays
- allocate(high_order_vert_flux(nVertLevels+1, nCells))
-
+ !$OMP DO PRIVATE(k)
do iCell = 1, nCells
do k=1, maxLevelCell(iCell)
inv_h_new(k, iCell) = 1.0 / (h(k, iCell) + dt * tend_h(k, iCell))
end do
end do
+ !$OMP END DO
! Loop over tracers. One tracer is advected at a time. It is copied into a temporary array in order to improve locality
+
do iTracer = 1, num_tracers
! Initialize variables for use in this iTracer iteration
+ !$OMP DO PRIVATE(k)
do iCell = 1, nCells
do k=1, maxLevelCell(iCell)
tracer_cur(k,iCell) = tracers(iTracer,k,iCell)
@@ -152,11 +161,15 @@
end if
end do ! k loop
end do ! iCell loop
+ !$OMP END DO
+ !$OMP WORKSHARE
high_order_vert_flux = 0.0
high_order_horiz_flux = 0.0
+ !$OMP END WORKSHARE
! Compute the high order vertical flux. Also determine bounds on tracer_cur.
+ !$OMP DO PRIVATE(k, i)
do iCell = 1, nCells
k = 1
tracer_max(k,iCell) = max(tracer_cur(k,iCell),tracer_cur(k+1,iCell))
@@ -192,8 +205,10 @@
end do ! k loop
end do ! i loop over nEdgesOnCell
end do ! iCell Loop
+ !$OMP END DO
! Compute the high order horizontal flux
+ !$OMP DO PRIVATE(i, iCell, k, tracer_weight)
do iEdge = 1, nEdges
do i = 1, nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
@@ -206,11 +221,13 @@
end do ! k loop
end do ! i loop over nAdvCellsForEdge
end do ! iEdge loop
+ !$OMP END DO
! low order upwind vertical flux (monotonic and diffused)
! Remove low order flux from the high order flux.
! Store left over high order flux in high_order_vert_flux array.
! Upwind fluxes are accumulated in upwind_tendency
+ !$OMP DO PRIVATE(k, flux_upwind)
do iCell = 1, nCells
do k = 2, maxLevelCell(iCell)
! dwj 02/03/12 Ocean and Atmosphere are different in vertical
@@ -234,11 +251,13 @@
flux_outgoing(k, iCell) = min(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - max(0.0_RKIND, high_order_vert_flux(k, iCell))
end do ! k Loop
end do ! iCell Loop
+ !$OMP END DO
! low order upwind horizontal flux (monotinc and diffused)
! Remove low order flux from the high order flux
! Store left over high order flux in high_order_horiz_flux array
! Upwind fluxes are accumulated in upwind_tendency
+ !$OMP DO PRIVATE(cell1, cell2, invAreaCell1, invAreaCell2, k, flux_upwind)
do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -251,7 +270,9 @@
high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) - flux_upwind
end do ! k loop
end do ! iEdge loop
+ !$OMP END DO
+ !$OMP DO PRIVATE(invAreaCell1, i, iEdge, cell1, cell2, k, flux_upwind)
do iCell = 1, nCells
invAreaCell1 = 1.0 / areaCell(iCell)
do i = 1, nEdgesOnCell(iCell)
@@ -269,10 +290,12 @@
end do
end do
end do
+ !$OMP END DO
! Build the factors for the FCT
! Computed using the bounds that were computed previously, and the bounds on the newly updated value
! Factors are placed in the flux_incoming and flux_outgoing arrays
+ !$OMP DO PRIVATE(k, tracer_min_new, tracer_max_new, tracer_upwind_new, scale_factor)
do iCell = 1, nCells
do k = 1, maxLevelCell(iCell)
tracer_min_new = (tracer_cur(k,iCell)*h(k,iCell) + dt*(upwind_tendency(k,iCell)+flux_outgoing(k,iCell))) * inv_h_new(k,iCell)
@@ -286,8 +309,10 @@
flux_outgoing(k,iCell) = min( 1.0_RKIND, max( 0.0_RKIND, scale_factor) )
end do ! k loop
end do ! iCell loop
+ !$OMP END DO
! rescale the high order horizontal fluxes
+ !$OMP DO PRIVATE(cell1, cell2, k, flux)
do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -298,8 +323,10 @@
high_order_horiz_flux(k,iEdge) = flux
end do ! k loop
end do ! iEdge loop
+ !$OMP END DO
! rescale the high order vertical flux
+ !$OMP DO PRIVATE(k, flux)
do iCell = 1, nCellsSolve
do k = 2, maxLevelCell(iCell)
flux = high_order_vert_flux(k,iCell)
@@ -311,8 +338,10 @@
high_order_vert_flux(k,iCell) = flux
end do ! k loop
end do ! iCell loop
+ !$OMP END DO
! Accumulate the scaled high order horizontal tendencies
+ !$OMP DO PRIVATE(invAreaCell1, i, iEdge, k)
do iCell = 1, nCells
invAreaCell1 = 1.0 / areaCell(iCell)
do i = 1, nEdgesOnCell(iCell)
@@ -326,8 +355,10 @@
end do
end do
end do
+ !$OMP END DO
! Accumulate the scaled high order vertical tendencies, and the upwind tendencies
+ !$OMP DO PRIVATE(k)
do iCell = 1, nCellsSolve
do k = 1,maxLevelCell(iCell)
tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
@@ -341,9 +372,11 @@
end if
end do ! k loop
end do ! iCell loop
+ !$OMP END DO
if (config_check_monotonicity) then
!build min and max bounds on old and new tracer for check on monotonicity.
+ !$OMP DO PRIVATE(k)
do iCell = 1, nCellsSolve
do k = 1, maxLevelCell(iCell)
if(tracer_new(k,iCell) < tracer_min(k, iCell)-eps) then
@@ -355,20 +388,34 @@
end if
end do
end do
+ !$OMP END DO
end if
end do ! iTracer loop
- deallocate(tracer_new)
- deallocate(tracer_cur)
- deallocate(upwind_tendency)
- deallocate(inv_h_new)
- deallocate(tracer_max)
- deallocate(tracer_min)
- deallocate(flux_incoming)
- deallocate(flux_outgoing)
- deallocate(high_order_horiz_flux)
- deallocate(high_order_vert_flux)
- deallocate(high_order_horiz_flux_field)
+
+ !$OMP SECTIONS
+ !$OMP SECTION
+ call mpas_deallocate_scratch_field(scratch % tracer_new)
+ !$OMP SECTION
+ call mpas_deallocate_scratch_field(scratch % tracer_cur)
+ !$OMP SECTION
+ call mpas_deallocate_scratch_field(scratch % upwind_tendency)
+ !$OMP SECTION
+ call mpas_deallocate_scratch_field(scratch % inv_h_new)
+ !$OMP SECTION
+ call mpas_deallocate_scratch_field(scratch % tracer_max)
+ !$OMP SECTION
+ call mpas_deallocate_scratch_field(scratch % tracer_min)
+ !$OMP SECTION
+ call mpas_deallocate_scratch_field(scratch % flux_incoming)
+ !$OMP SECTION
+ call mpas_deallocate_scratch_field(scratch % flux_outgoing)
+ !$OMP SECTION
+ call mpas_deallocate_scratch_field(scratch % high_order_horiz_flux)
+ !$OMP SECTION
+ call mpas_deallocate_scratch_field(scratch % high_order_vert_flux)
+ !$OMP END SECTIONS
+
end subroutine mpas_ocn_tracer_advection_mono_tend!}}}
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -61,12 +61,18 @@
real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !< Input: Distance between vertical interfaces of a cell
type (mesh_type), intent(in) :: grid !< Input: Grid information
+ !$OMP SINGLE
call mpas_timer_start("tracer-hadv", .false.)
+ !$OMP END SINGLE
call mpas_ocn_tracer_advection_std_hadv_tend(tracers, uh, grid, tend)
+ !$OMP SINGLE
call mpas_timer_stop("tracer-hadv")
call mpas_timer_start("tracer-vadv", .false.)
+ !$OMP END SINGLE
call mpas_ocn_tracer_advection_std_vadv_tend(tracers, w, verticalCellSize, grid, tend)
+ !$OMP SINGLE
call mpas_timer_stop("tracer-vadv")
+ !$OMP END SINGLE
end subroutine mpas_ocn_tracer_advection_std_tend!}}}
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std_hadv.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std_hadv.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std_hadv.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -82,6 +82,7 @@
allocate(flux_arr(num_tracers, grid % nVertLevels))
! horizontal flux divergence, accumulate in tracer_tend
+ !$OMP DO PRIVATE(cell1, cell2, i, iCell, k, tracer_weight, iTracer)
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -106,6 +107,7 @@
end do
end if
end do
+ !$OMP END DO
deallocate(flux_arr)
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std_vadv2.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std_vadv2.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std_vadv2.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -71,6 +71,7 @@
vert_flux(:,1) = 0.
+ !$OMP DO PRIVATE(k, iTracer, weightK, weightKm1)
do iCell=1,grid % nCellsSolve
do k = 2, maxLevelCell(iCell)
do iTracer=1,num_tracers
@@ -88,6 +89,7 @@
end do
end do
end do
+ !$OMP END DO
deallocate(vert_flux)
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std_vadv3.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std_vadv3.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std_vadv3.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -65,6 +65,7 @@
vert_flux(:,1) = 0.
+ !$OMP DO PRIVATE(k, iTracer, weightK, weightKm1)
do iCell=1,grid % nCellsSolve
k = 2
@@ -98,6 +99,7 @@
end do
end do
end do
+ !$OMP END DO
deallocate(vert_flux)
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std_vadv4.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std_vadv4.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_advection_std_vadv4.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -66,6 +66,7 @@
vert_flux(:,1) = 0.
+ !$OMP DO PRIVATE(k, iTracer, weightK, weightKm1)
do iCell=1,grid % nCellsSolve
k = 2
@@ -95,8 +96,8 @@
tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (vert_flux(iTracer, k+1) - vert_flux(iTracer, k))
end do
end do
-
end do
+ !$OMP END DO
deallocate(vert_flux)
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_hmix.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_hmix.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_hmix.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -72,7 +72,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tracer_hmix_tend(grid, h_edge, tracers, tend, err)!{{{
+ subroutine ocn_tracer_hmix_tend(grid, h_edge, tracers, tend, scratch, err)!{{{
!-----------------------------------------------------------------
!
@@ -98,6 +98,8 @@
real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
tend !< Input/Output: velocity tendency
+ type (scratch_type), intent(inout) :: scratch !< Input/Output: Scratch structure
+
!-----------------------------------------------------------------
!
! output variables
@@ -122,12 +124,18 @@
!
!-----------------------------------------------------------------
+ !$OMP SINGLE
call mpas_timer_start("del2", .false., del2Timer)
+ !$OMP END SINGLE
call ocn_tracer_hmix_del2_tend(grid, h_edge, tracers, tend, err1)
+ !$OMP SINGLE
call mpas_timer_stop("del2", del2Timer)
call mpas_timer_start("del4", .false., del4Timer)
- call ocn_tracer_hmix_del4_tend(grid, h_edge, tracers, tend, err2)
+ !$OMP END SINGLE
+ call ocn_tracer_hmix_del4_tend(grid, h_edge, tracers, tend, scratch, err2)
+ !$OMP SINGLE
call mpas_timer_stop("del4", del4Timer)
+ !$OMP END SINGLE
err = ior(err1, err2)
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_hmix_del2.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -153,6 +153,7 @@
!
! compute a boundary mask to enforce insulating boundary conditions in the horizontal
!
+ !$OMP DO PRIVATE(invAreaCell1, i, iEdge, cell1, cell2, r_tmp, k, iTracer, flux)
do iCell = 1, nCells
invAreaCell1 = 1.0 / areaCell(iCell)
do i = 1, nEdgesOncell(iCell)
@@ -176,6 +177,7 @@
end do
end do
+ !$OMP END DO
!--------------------------------------------------------------------
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_hmix_del4.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -67,7 +67,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tracer_hmix_del4_tend(grid, h_edge, tracers, tend, err)!{{{
+ subroutine ocn_tracer_hmix_del4_tend(grid, h_edge, tracers, tend, scratch, err)!{{{
!-----------------------------------------------------------------
!
@@ -93,6 +93,8 @@
real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
tend !< Input/Output: velocity tendency
+ type (scratch_type), intent(inout) :: scratch !< Scratch structure
+
!-----------------------------------------------------------------
!
! output variables
@@ -115,7 +117,7 @@
real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux, flux, invdcEdge, r_tmp1, r_tmp2
- real (kind=RKIND), dimension(:,:,:), allocatable :: delsq_tracer
+ real (kind=RKIND), dimension(:,:,:), pointer :: delsq_tracer
real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaCell, meshScalingDel4
@@ -152,11 +154,19 @@
edgesOnCell => grid % edgesOnCell % array
edgeSignOnCell => grid % edgeSignOnCell % array
- allocate(delsq_tracer(num_tracers,nVertLevels, nCells+1))
+ !$OMP SINGLE
+ call mpas_allocate_scratch_field(scratch % delsq_tracers)
+! allocate(delsq_tracer(num_tracers,nVertLevels, nCells+1))
+ !$OMP END SINGLE
+
+ delsq_tracer => scratch % delsq_tracers % array
+ !$OMP WORKSHARE
delsq_tracer(:,:,:) = 0.0
+ !$OMP END WORKSHARE
! first del2: div(h </font>
<font color="gray">abla \phi) at cell center
+ !$OMP DO PRIVATE(invAreaCell1, i, iEdge, invdcEdge, cell1, cell2, k, iTracer, r_tmp1, r_tmp2)
do iCell = 1, nCells
invAreaCell1 = 1.0 / areaCell(iCell)
do i = 1, nEdgesOnCell(iCell)
@@ -176,8 +186,10 @@
end do
end do
end do
+ !$OMP END DO
! second del2: div(h </font>
<font color="gray">abla [delsq_tracer]) at cell center
+ !$OMP DO PRIVATE(invAreaCell1, i, iEdge, cell1, cell2, invdcEdge, k, iTracer, tracer_turb_flux, flux)
do iCell = 1, nCells
invAreaCell1 = 1.0 / areaCell(iCell)
do i = 1, nEdgesOnCell(iCell)
@@ -198,8 +210,12 @@
end do
end do
end do
+ !$OMP END DO
- deallocate(delsq_tracer)
+ !$OMP SINGLE
+ call mpas_deallocate_scratch_field(scratch % delsq_tracers)
+! deallocate(delsq_tracer)
+ !$OMP END SINGLE
!--------------------------------------------------------------------
end subroutine ocn_tracer_hmix_del4_tend!}}}
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_coriolis.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_coriolis.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_coriolis.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -127,6 +127,7 @@
nEdgesSolve = grid % nEdgesSolve
+ !$OMP DO PRIVATE(cell1, cell2, invLength, k, q, j, eoe, workpv)
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -146,6 +147,7 @@
end do
end do
+ !$OMP END DO
!--------------------------------------------------------------------
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -127,6 +127,7 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
edgeMask => grid % edgeMask % array
+ !$OMP DO PRIVATE(k)
do iEdge=1,grid % nEdgesSolve
k = max(maxLevelEdgeTop(iEdge), 1)
@@ -138,6 +139,7 @@
tend(k,iEdge) = tend(k,iEdge)-edgeMask(k,iEdge)*(bottomDragCoef*u(k,iEdge)*sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge))
enddo
+ !$OMP END DO
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_forcing_rayleigh.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_forcing_rayleigh.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_forcing_rayleigh.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -119,6 +119,7 @@
nEdgesSolve = grid % nEdgesSolve
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ !$OMP DO PRIVATE(k)
do iEdge=1,nEdgesSolve
do k=1,maxLevelEdgeTop(iEdge)
@@ -126,8 +127,8 @@
enddo
enddo
+ !$OMP END DO
-
!--------------------------------------------------------------------
end subroutine ocn_vel_forcing_rayleigh_tend!}}}
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_forcing_windstress.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_forcing_windstress.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_forcing_windstress.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -122,6 +122,7 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
edgeMask => grid % edgeMask % array
+ !$OMP DO PRIVATE(k)
do iEdge=1,nEdgesSolve
! efficiency note: it would be nice to avoid this
! if within a do. This could be done with
@@ -132,8 +133,8 @@
! forcing in top layer only
tend(k,iEdge) = tend(k,iEdge) + edgeMask(k, iEdge) * (u_src(k,iEdge) / config_rho0 / h_edge(k,iEdge))
enddo
-
enddo
+ !$OMP END DO
!--------------------------------------------------------------------
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_hmix.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_hmix.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_hmix.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -73,7 +73,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_vel_hmix_tend(grid, divergence, vorticity, viscosity, tend, err)!{{{
+ subroutine ocn_vel_hmix_tend(grid, divergence, vorticity, viscosity, tend, scratch, err)!{{{
!-----------------------------------------------------------------
!
@@ -100,8 +100,10 @@
tend !< Input/Output: velocity tendency
real (kind=RKIND), dimension(:,:), intent(inout) :: &
- viscosity !< Input: viscosity
+ viscosity !< Input/Output: viscosity
+ type (scratch_type), intent(inout) :: scratch !< Input/Output: Scratch structure
+
!-----------------------------------------------------------------
!
! output variables
@@ -137,7 +139,7 @@
call mpas_timer_stop("leith", leithTimer)
call mpas_timer_start("del4", .false., del4Timer)
- call ocn_vel_hmix_del4_tend(grid, divergence, vorticity, tend, err3)
+ call ocn_vel_hmix_del4_tend(grid, divergence, vorticity, tend, scratch, err3)
call mpas_timer_stop("del4", del4Timer)
err = ior(ior(err1, err2),err3)
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_hmix_del2.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_hmix_del2.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -140,6 +140,7 @@
dcEdge => grid % dcEdge % array
dvEdge => grid % dvEdge % array
+ !$OMP DO PRIVATE(cell1, cell2, vertex1, vertex2, invLength1, invLength2, k, u_diffusion, visc2)
do iEdge=1,nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -167,6 +168,7 @@
end do
end do
+ !$OMP END DO
!--------------------------------------------------------------------
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_hmix_del4.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -72,7 +72,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_vel_hmix_del4_tend(grid, divergence, vorticity, tend, err)!{{{
+ subroutine ocn_vel_hmix_del4_tend(grid, divergence, vorticity, tend, scratch, err)!{{{
!-----------------------------------------------------------------
!
@@ -98,6 +98,8 @@
real (kind=RKIND), dimension(:,:), intent(inout) :: &
tend !< Input/Output: velocity tendency
+ type (scratch_type), intent(inout) :: scratch !< Input/Output: Scratch structure
+
!-----------------------------------------------------------------
!
! output variables
@@ -126,8 +128,8 @@
real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaTriangle, &
meshScalingDel4, areaCell
- real (kind=RKIND), dimension(:,:), allocatable :: delsq_divergence, &
- delsq_circulation, delsq_vorticity, delsq_u
+ real (kind=RKIND), dimension(:,:), pointer :: delsq_divergence, &
+ delsq_vorticity, delsq_u
err = 0
@@ -156,16 +158,28 @@
edgesOnCell => grid % edgesOnCell % array
edgeSignOnVertex => grid % edgeSignOnVertex % array
edgeSignOnCell => grid % edgeSignOnCell % array
+ !DWJ 01/29/13 OMP DEL4
+ !$OMP MASTER
+ call mpas_allocate_scratch_field(scratch % delsq_u, .true.)
+ call mpas_allocate_scratch_field(scratch % delsq_divergence, .true.)
+ call mpas_allocate_scratch_field(scratch % delsq_vorticity, .true.)
+ !$OMP END MASTER
- allocate(delsq_u(nVertLEvels, nEdges+1))
- allocate(delsq_divergence(nVertLevels, nCells+1))
- allocate(delsq_vorticity(nVertLevels, nVertices+1))
+ !$OMP BARRIER
+ delsq_u => scratch % delsq_u % array
+ delsq_divergence => scratch % delsq_divergence % array
+ delsq_vorticity => scratch % delsq_vorticity % array
+
+
+ !$OMP WORKSHARE
delsq_u(:,:) = 0.0
delsq_vorticity(:,:) = 0.0
delsq_divergence(:,:) = 0.0
+ !$OMP END WORKSHARE
!Compute delsq_u
+ !$OMP DO PRIVATE(cell1, cell2, vertex1, vertex2, invDcEdge, invDvEdge, k)
do iEdge = 1, nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -182,8 +196,10 @@
-viscVortCoef *( vorticity(k,vertex2) - vorticity(k,vertex1)) * invDcEdge * sqrt(3.0) ! TDR
end do
end do
+ !$OMP END DO
! Compute delsq_vorticity
+ !$OMP DO PRIVATE(invAreaTri1, i, iEdge, k)
do iVertex = 1, nVertices
invAreaTri1 = 1.0 / areaTriangle(iVertex)
do i = 1, vertexDegree
@@ -193,8 +209,10 @@
end do
end do
end do
+ !$OMP END DO
! Compute delsq_divergence
+ !$OMP DO PRIVATE(invAreaCell1, i, iEdge, k)
do iCell = 1, nCells
invAreaCell1 = 1.0 / areaCell(iCell)
do i = 1, nEdgesOnCell(iCell)
@@ -204,9 +222,11 @@
end do
end do
end do
+ !$OMP END DO
! Compute - \kappa </font>
<font color="black">abla^4 u
! as </font>
<font color="black">abla div(</font>
<font color="black">abla^2 u) + k \times </font>
<font color="black">abla ( k \cross curl(</font>
<font color="gray">abla^2 u) )
+ !$OMP DO PRIVATE(cell1, cell2, vertex1, vertex2, invDcEdge, invDvEdge, r_tmp, k, u_diffusion)
do iEdge=1,nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -224,10 +244,14 @@
tend(k,iEdge) = tend(k,iEdge) - edgeMask(k, iEdge) * u_diffusion * r_tmp
end do
end do
+ !$OMP END DO
- deallocate(delsq_u)
- deallocate(delsq_divergence)
- deallocate(delsq_vorticity)
+ !$OMP BARRIER
+ !$OMP SINGLE
+ call mpas_deallocate_scratch_field(scratch % delsq_u, .true.)
+ call mpas_deallocate_scratch_field(scratch % delsq_divergence, .true.)
+ call mpas_deallocate_scratch_field(scratch % delsq_vorticity, .true.)
+ !$OMP END SINGLE
!--------------------------------------------------------------------
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_pressure_grad.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_pressure_grad.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_pressure_grad.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -127,6 +127,7 @@
! we have set rho0Inv=1 and grho0Inv=0 in the init routine,
! and pressure is passed in as MontPot.
+ !$OMP DO PRIVATE(cell1, cell2, invdcEdge, k)
do iEdge=1,nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -141,8 +142,8 @@
- zMid(k,cell1) )* invdcEdge
end do
-
end do
+ !$OMP END DO
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_vadv.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_vadv.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vel_vadv.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -129,6 +129,7 @@
allocate(w_dudzTopEdge(nVertLevels+1))
w_dudzTopEdge = 0.0
+ !$OMP DO PRIVATE(cell1, cell2, k, wTopEdge)
do iEdge=1,nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -148,6 +149,7 @@
tend(k,iEdge) = tend(k,iEdge) - edgeMask(k, iEdge) * 0.5 * (w_dudzTopEdge(k) + w_dudzTopEdge(k+1))
enddo
enddo
+ !$OMP END DO
deallocate(w_dudzTopEdge)
!--------------------------------------------------------------------
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vmix.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vmix.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vmix.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -76,7 +76,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_vmix_coefs(grid, s, d, err)!{{{
+ subroutine ocn_vmix_coefs(grid, s, d, scratch, err)!{{{
!-----------------------------------------------------------------
!
@@ -99,6 +99,8 @@
type (diagnostics_type), intent(inout) :: &
d !< Input/Output: diagnostic information
+ type (scratch_type), intent(inout) :: scratch !< Input/Output: scratch structure
+
!-----------------------------------------------------------------
!
! output variables
@@ -125,7 +127,7 @@
call ocn_vmix_coefs_const_build(grid, s, d, err1)
call ocn_vmix_coefs_tanh_build(grid, s, d, err2)
- call ocn_vmix_coefs_rich_build(grid, s, d, err3)
+ call ocn_vmix_coefs_rich_build(grid, s, d, scratch, err3)
err = ior(err1, ior(err2, err3))
@@ -309,6 +311,7 @@
allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),uTemp(nVertLevels))
A(1)=0
+ !$OMP DO PRIVATE(n, cell1, cell2, k)
do iEdge=1,nEdges
N=maxLevelEdgeTop(iEdge)
if (N.gt.0) then
@@ -353,6 +356,7 @@
end if
end do
+ !$OMP END DO
deallocate(A,B,C,uTemp)
@@ -539,6 +543,7 @@
allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),tracersTemp(num_tracers,nVertLevels))
+ !$OMP DO PRIVATE(n, k)
do iCell=1,nCells
! Compute A(k), B(k), C(k) for tracers
N = maxLevelCell(iCell)
@@ -568,6 +573,7 @@
tracers(:,1:N,iCell) = tracersTemp(:,1:N)
tracers(:,N+1:nVertLevels,iCell) = -1e34
end do
+ !$OMP END DO
deallocate(A,B,C,tracersTemp)
@@ -590,11 +596,12 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_vmix_implicit(dt, grid, diagnostics, state, err)!{{{
+ subroutine ocn_vmix_implicit(dt, grid, diagnostics, state, scratch, err)!{{{
real (kind=RKIND), intent(in) :: dt
type (mesh_type), intent(in) :: grid
type (diagnostics_type), intent(inout) :: diagnostics
type (state_type), intent(inout) :: state
+ type (scratch_Type), intent(inout) :: scratch
integer, intent(out) :: err
integer :: nCells
@@ -615,7 +622,7 @@
nCells = grid % nCells
- call ocn_vmix_coefs(grid, state, diagnostics, err)
+ call ocn_vmix_coefs(grid, state, diagnostics, scratch, err)
!
! Implicit vertical solve for momentum
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vmix_coefs_const.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vmix_coefs_const.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vmix_coefs_const.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -185,7 +185,9 @@
if(.not.constViscOn) return
+ !$OMP WORKSHARE
vertViscTopOfEdge = constVisc
+ !$OMP END WORKSHARE
!--------------------------------------------------------------------
@@ -241,7 +243,9 @@
if(.not.constDiffOn) return
+ !$OMP WORKSHARE
vertDiffTopOfCell = constDiff
+ !$OMP END WORKSHARE
!--------------------------------------------------------------------
Modified: branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2013-01-31 17:05:58 UTC (rev 2403)
@@ -69,7 +69,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_vmix_coefs_rich_build(grid, s, d, err)!{{{
+ subroutine ocn_vmix_coefs_rich_build(grid, s, d, scratch, err)!{{{
!-----------------------------------------------------------------
!
@@ -92,6 +92,8 @@
type (diagnostics_type), intent(inout) :: &
d !< Input/Output: diagnostic information
+ type (scratch_type), intent(inout) :: scratch !< Input/Output: Scratch structure
+
!-----------------------------------------------------------------
!
! output variables
@@ -141,13 +143,15 @@
rhoDisplaced => s % rhoDisplaced % array
tracers => s % tracers % array
+ !$OMP SINGLE
call mpas_timer_start("eos rich", .false., richEOSTimer)
call ocn_equation_of_state_rho(s, grid, 0,'relative', err)
call ocn_equation_of_state_rho(s, grid, 1,'relative', err)
call mpas_timer_stop("eos rich", richEOSTimer)
+ !$OMP END SINGLE
call ocn_vmix_get_rich_numbers(grid, indexT, indexS, u, h, h_edge, &
- rho, rhoDisplaced, tracers, RiTopOfEdge, RiTopOfCell, err1)
+ rho, rhoDisplaced, tracers, RiTopOfEdge, RiTopOfCell, scratch, err1)
call ocn_vel_vmix_coefs_rich(grid, RiTopOfEdge, h_edge, vertViscTopOfEdge, err2)
call ocn_tracer_vmix_coefs_rich(grid, RiTopOfCell, h, vertDiffTopOfCell, err3)
@@ -223,6 +227,7 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
vertViscTopOfEdge = 0.0
+ !$OMP DO PRIVATE(k)
do iEdge = 1,nEdges
do k = 2,maxLevelEdgeTop(iEdge)
! mrp 110324 efficiency note: this if is inside iEdge and k loops.
@@ -254,6 +259,7 @@
end if
end do
end do
+ !$OMP END DO
!--------------------------------------------------------------------
@@ -328,6 +334,7 @@
vertDiffTopOfCell = 0.0
coef = -gravity/config_rho0/2.0
+ !$OMP DO PRIVATE(k)
do iCell = 1,nCells
do k = 2,maxLevelCell(iCell)
! mrp 110324 efficiency note: this if is inside iCell and k loops.
@@ -361,8 +368,8 @@
end if
end do
end do
+ !$OMP END DO
-
!--------------------------------------------------------------------
end subroutine ocn_tracer_vmix_coefs_rich!}}}
@@ -382,7 +389,7 @@
!-----------------------------------------------------------------------
subroutine ocn_vmix_get_rich_numbers(grid, indexT, indexS, u, h, h_edge, & !{{{
- rho, rhoDisplaced, tracers, RiTopOfEdge, RiTopOfCell, err)
+ rho, rhoDisplaced, tracers, RiTopOfEdge, RiTopOfCell, scratch, err)
!-----------------------------------------------------------------
!
@@ -419,6 +426,8 @@
real (kind=RKIND), dimension(:,:), intent(inout) :: RiTopOfEdge !< Input/output: Richardson number top of cell
real (kind=RKIND), dimension(:,:), intent(inout) :: RiTopOfCell !< Input/output: Richardson number top of cell
+ type (scratch_type), intent(inout) :: scratch !< Input/Output: Scratch structure
+
integer, intent(inout) :: err !< Output: error flag
!-----------------------------------------------------------------
@@ -435,7 +444,7 @@
real (kind=RKIND) :: coef, invAreaCell
real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaCell
- real (kind=RKIND), dimension(:,:), allocatable :: drhoTopOfCell, du2TopOfCell, &
+ real (kind=RKIND), dimension(:,:), pointer :: drhoTopOfCell, du2TopOfCell, &
drhoTopOfEdge, du2TopOfEdge
err = 0
@@ -457,9 +466,19 @@
edgesOnCell => grid % edgesOnCell % array
edgeSignOnCell => grid % edgeSignOnCell % array
- allocate( &
- drhoTopOfCell(nVertLevels+1,nCells+1), drhoTopOfEdge(nVertLevels+1,nEdges), &
- du2TopOfCell(nVertLevels+1,nCells+1), du2TopOfEdge(nVertLevels+1,nEdges))
+ !$OMP SINGLE
+ call mpas_allocate_scratch_field(scratch % drhoTopOfCell)
+ call mpas_allocate_scratch_field(scratch % drhoTopOfEdge)
+ call mpas_allocate_scratch_field(scratch % du2TopOfCell)
+ call mpas_allocate_scratch_field(scratch % du2TopOfEdge)
+! allocate( &
+! drhoTopOfCell(nVertLevels+1,nCells+1), drhoTopOfEdge(nVertLevels+1,nEdges), &
+! du2TopOfCell(nVertLevels+1,nCells+1), du2TopOfEdge(nVertLevels+1,nEdges))
+ !$OMP END SINGLE
+ drhoTopOfCell => scratch % drhoTopOfCell % array
+ drhoTopOfEdge => scratch % drhoTopOfEdge % array
+ du2TopOfCell => scratch % du2TopOfCell % array
+ du2TopOfEdge => scratch % du2TopOfEdge % array
! compute density of parcel displaced to next deeper z-level,
! in state % rhoDisplaced
@@ -472,13 +491,19 @@
! drhoTopOfCell(k) = $\rho^*_{k-1}-\rho^*_k$
+ !$OMP WORKSHARE
drhoTopOfCell = 0.0
+ !$OMP END WORKSHARE
+
+ !$OMP DO PRIVATE(k)
do iCell=1,nCells
do k=2,maxLevelCell(iCell)
drhoTopOfCell(k,iCell) = rhoDisplaced(k-1,iCell) - rhoDisplaced(k,iCell)
end do
end do
+ !$OMP END DO
+ !$OMP SINGLE
! interpolate drhoTopOfCell to drhoTopOfEdge
drhoTopOfEdge = 0.0
do iEdge=1,nEdges
@@ -534,9 +559,16 @@
/ (du2TopOfCell(k,iCell) + 1e-20)
end do
end do
+ !$OMP END SINGLE
- deallocate(drhoTopOfCell, drhoTopOfEdge, &
- du2TopOfCell, du2TopOfEdge)
+ !$OMP SINGLE
+ call mpas_allocate_scratch_field(scratch % drhoTopOfCell)
+ call mpas_allocate_scratch_field(scratch % drhoTopOfEdge)
+ call mpas_allocate_scratch_field(scratch % du2TopOfCell)
+ call mpas_allocate_scratch_field(scratch % du2TopOfEdge)
+! deallocate(drhoTopOfCell, drhoTopOfEdge, &
+! du2TopOfCell, du2TopOfEdge)
+ !$OMP END SINGLE
!--------------------------------------------------------------------
Modified: branches/ocean_projects/openmp_elements/src/registry/gen_inc.c
===================================================================
--- branches/ocean_projects/openmp_elements/src/registry/gen_inc.c 2013-01-31 17:02:24 UTC (rev 2402)
+++ branches/ocean_projects/openmp_elements/src/registry/gen_inc.c 2013-01-31 17:05:58 UTC (rev 2403)
@@ -718,35 +718,37 @@
var_list_ptr3 = var_list_ptr3->next;
}
- fortprintf(fd, " allocate(%s %% %s %% array(%i, ", group_ptr->name, var_ptr2->super_array, i);
- dimlist_ptr = var_ptr2->dimlist;
- if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
- !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
- !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
- if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_code);
- else fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_file);
- else
- if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, "%s", dimlist_ptr->dim->name_in_file);
- else fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
- dimlist_ptr = dimlist_ptr->next;
- while (dimlist_ptr) {
+ if(var_ptr2->persistence == PERSISTENT){
+ fortprintf(fd, " allocate(%s %% %s %% array(%i, ", group_ptr->name, var_ptr2->super_array, i);
+ dimlist_ptr = var_ptr2->dimlist;
if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
!strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
- if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_code);
- else fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_file);
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, "%s + 1", dimlist_ptr->dim->name_in_file);
else
- if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_file);
- else fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, "%s", dimlist_ptr->dim->name_in_file);
+ else fortprintf(fd, "%s", dimlist_ptr->dim->name_in_code);
dimlist_ptr = dimlist_ptr->next;
- }
- fortprintf(fd, "))</font>
<font color="red">");
- if (var_ptr->vtype == INTEGER)
- fortprintf(fd, " %s %% %s %% array = 0</font>
<font color="red">", group_ptr->name, var_ptr2->super_array ); /* initialize field to zero */
- else if (var_ptr->vtype == REAL)
- fortprintf(fd, " %s %% %s %% array = 0.0</font>
<font color="red">", group_ptr->name, var_ptr2->super_array ); /* initialize field to zero */
- else if (var_ptr->vtype == CHARACTER)
- fortprintf(fd, " %s %% %s %% array = \'\'</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array ); /* initialize field to zero */
+ while (dimlist_ptr) {
+ if (!strncmp(dimlist_ptr->dim->name_in_file, "nCells", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
+ !strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
+ if (!dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_code);
+ else fortprintf(fd, ", %s + 1", dimlist_ptr->dim->name_in_file);
+ else
+ if (dimlist_ptr->dim->namelist_defined) fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_file);
+ else fortprintf(fd, ", %s", dimlist_ptr->dim->name_in_code);
+ dimlist_ptr = dimlist_ptr->next;
+ }
+ fortprintf(fd, "))</font>
<font color="blue">");
+ if (var_ptr->vtype == INTEGER)
+ fortprintf(fd, " %s %% %s %% array = 0</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array ); /* initialize field to zero */
+ else if (var_ptr->vtype == REAL)
+ fortprintf(fd, " %s %% %s %% array = 0.0</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array ); /* initialize field to zero */
+ else if (var_ptr->vtype == CHARACTER)
+ fortprintf(fd, " %s %% %s %% array = \'\'</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array ); /* initialize field to zero */
+ }
fortprintf(fd, " %s %% %s %% dimSizes(1) = %i</font>
<font color="black">", group_ptr->name, var_ptr2->super_array, i);
fortprintf(fd, " %s %% %s %% dimNames(1) = \'num_%s\'</font>
<font color="gray">", group_ptr->name, var_ptr2->super_array, var_ptr2->super_array);
@@ -757,8 +759,14 @@
!strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
if (!dimlist_ptr->dim->namelist_defined) {
- fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s</font>
<font color="red">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_code);
- fortprintf(fd, " %s %% %s %% dimNames(%i) = \'%s\'</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_file);
+ if (var_ptr2->persistence == PERSISTENT){
+ fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_code);
+ fortprintf(fd, " %s %% %s %% dimNames(%i) = \'%s\'</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_file);
+ }
+ else {
+ fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s+1</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_code);
+ fortprintf(fd, " %s %% %s %% dimNames(%i) = \'%s\'</font>
<font color="blue">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_file);
+ }
}
else {
fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s</font>
<font color="gray">", group_ptr->name, var_ptr2->super_array, i, dimlist_ptr->dim->name_in_file);
@@ -814,7 +822,6 @@
fortprintf(fd, " %s %% %s %% isSuperArray = .false.</font>
<font color="red">", group_ptr->name, var_ptr->name_in_code);
if (var_ptr->ndims > 0) {
if(var_ptr->persistence == SCRATCH){
- fortprintf(fd, " ! SCRATCH VARIABLE</font>
<font color="black">");
fortprintf(fd, " nullify(%s %% %s %% array)</font>
<font color="gray">", group_ptr->name, var_ptr->name_in_code);
} else if(var_ptr->persistence == PERSISTENT){
fortprintf(fd, " allocate(%s %% %s %% array(", group_ptr->name, var_ptr->name_in_code);
@@ -855,8 +862,14 @@
!strncmp(dimlist_ptr->dim->name_in_file, "nEdges", 1024) ||
!strncmp(dimlist_ptr->dim->name_in_file, "nVertices", 1024))
if (!dimlist_ptr->dim->namelist_defined) {
- fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s</font>
<font color="red">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_code);
- fortprintf(fd, " %s %% %s %% dimNames(%i) = \'%s\'</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_file);
+ if(var_ptr->persistence == PERSISTENT){
+ fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_code);
+ fortprintf(fd, " %s %% %s %% dimNames(%i) = \'%s\'</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_file);
+ }
+ else {
+ fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s+1</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_code);
+ fortprintf(fd, " %s %% %s %% dimNames(%i) = \'%s\'</font>
<font color="blue">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_file);
+ }
}
else {
fortprintf(fd, " %s %% %s %% dimSizes(%i) = %s</font>
<font color="gray">", group_ptr->name, var_ptr->name_in_code, i, dimlist_ptr->dim->name_in_file);
@@ -999,10 +1012,12 @@
fortprintf(fd, " dest %% %s %% scalar = src %% %s %% scalar</font>
<font color="blue">", var_ptr2->super_array, var_ptr2->super_array);
}
else {
+ if (var_ptr->persistence == PERSISTENT){
if (var_ptr->ndims > 0)
fortprintf(fd, " dest %% %s %% array = src %% %s %% array</font>
<font color="black">", var_ptr->name_in_code, var_ptr->name_in_code);
else
fortprintf(fd, " dest %% %s %% scalar = src %% %s %% scalar</font>
<font color="blue">", var_ptr->name_in_code, var_ptr->name_in_code);
+ }
var_list_ptr = var_list_ptr->next;
}
}
</font>
</pre>