<p><b>dwj07@fsu.edu</b> 2012-08-24 15:18:43 -0600 (Fri, 24 Aug 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
3rd (?) cut at openmp on element loops.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/omp_blocks/openmp_test/src/core_sw_elements/Registry
===================================================================
--- branches/omp_blocks/openmp_test/src/core_sw_elements/Registry 2012-08-23 21:56:35 UTC (rev 2122)
+++ branches/omp_blocks/openmp_test/src/core_sw_elements/Registry 2012-08-24 21:18:43 UTC (rev 2123)
@@ -164,3 +164,14 @@
var persistent real gradPVn ( nVertLevels nEdges Time ) 2 - gradPVn state - -
var persistent real h_vertex ( nVertLevels nVertices Time ) 2 - h_vertex state - -
+% Test global arrays
+var persistent real delsq_divergence ( nVertLevels nCells ) 2 - delsq_divergence state - -
+var persistent real delsq_u ( nVertLevels nEdges ) 2 - delsq_u state - -
+var persistent real delsq_circulation ( nVertLevels nVertices ) 2 - delsq_circulation state - -
+var persistent real delsq_vorticity ( nVertLevels nVertices ) 2 - delsq_vorticity state - -
+var persistent real boundaryMask ( nVertLevels nEdges ) 0 - boundaryMask mesh - -
+var persistent real delsq_tracer ( nTracers nVertLevels nCells ) 2 - delsq_tracer state - -
+
+% Test Sign Arrays
+var persistent integer edgeSignOnCell ( maxEdges nCells ) 0 - edgeSignOnCell mesh - -
+var persistent integer edgeSignOnVertex ( maxEdges nVertices ) 0 - edgeSignOnVertex mesh - -
Modified: branches/omp_blocks/openmp_test/src/core_sw_elements/mpas_sw_mpas_core.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_sw_elements/mpas_sw_mpas_core.F 2012-08-23 21:56:35 UTC (rev 2122)
+++ branches/omp_blocks/openmp_test/src/core_sw_elements/mpas_sw_mpas_core.F 2012-08-24 21:18:43 UTC (rev 2123)
@@ -125,6 +125,7 @@
call sw_compute_solve_diagnostics(dt, block % state % time_levs(1) % state, mesh)
call compute_mesh_scaling(mesh)
+ call setup_sign_arrays(mesh)
call mpas_rbf_interp_initialize(mesh)
call mpas_init_reconstruct(mesh)
@@ -163,8 +164,6 @@
integer :: ierr
! Eventually, dt should be domain specific
- !$omp parallel default(shared)
- !$omp single
dt = config_dt
currTime = mpas_get_clock_time(clock, MPAS_NOW, ierr)
@@ -176,10 +175,8 @@
! 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
- !$omp end single
do while (.not. mpas_is_clock_stop_time(clock))
- !$omp single
itimestep = itimestep + 1
call mpas_advance_clock(clock)
@@ -188,28 +185,19 @@
write(0,*) 'Doing timestep ', trim(timeStamp)
call mpas_timer_start("time integration")
- !$omp end single
call mpas_timestep(domain, itimestep, dt, timeStamp)
- !$omp single
call mpas_timer_stop("time integration")
- !$omp end single
! Move time level 2 fields back into time level 1 for next time step
- !$omp single private(block_ptr)
block_ptr => domain % blocklist
do while(associated(block_ptr))
block = block_ptr
- !$omp task firstprivate(block)
call mpas_shift_time_levels_state(block % state)
- !$omp end task
block_ptr => block_ptr % next
end do
- !$omp taskwait
- !$omp end single
!TODO: mpas_get_clock_ringing_alarms is probably faster than multiple mpas_is_alarm_ringing...
- !$omp single
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 maximum number of frames per outfile was reached
@@ -228,10 +216,8 @@
call mpas_output_state_for_domain(restart_obj, domain, 1)
call mpas_output_state_finalize(restart_obj, domain % dminfo)
end if
- !$omp end single
end do
- !$omp end parallel
end subroutine mpas_core_run
@@ -401,4 +387,40 @@
end subroutine compute_mesh_scaling
+ subroutine setup_sign_arrays(mesh)
+
+ use mpas_grid_types
+
+ implicit none
+
+ type (mesh_type), intent(inout) :: mesh
+
+ integer :: iCell, iVertex, iEdge, i
+
+ do iCell = 1, mesh % nCells
+ do i = 1, mesh % nEdgesOnCell % array(iCell)
+ iEdge = mesh % edgesOnCell % array(i, iCell)
+
+ if(mesh % cellsOnEdge % array(1, iEdge) == iCell) then
+ mesh % edgeSignOnCell % array(i, iCell) = -1
+ else
+ mesh % edgeSignOnCell % array(i, iCell) = 1
+ end if
+ end do
+ end do
+
+ do iVertex = 1, mesh % nVertices
+ do i = 1, mesh % vertexDegree
+ iEdge = mesh % edgesOnVertex % array(i, iVertex)
+
+ if(mesh % verticesOnEdge % array(1, iEdge) == iVertex) then
+ mesh % edgeSignOnVertex % array(i, iVertex) = -1
+ else
+ mesh % edgeSignOnVertex % array(i, iVertex) = 1
+ end if
+ end do
+ end do
+
+ end subroutine setup_sign_arrays
+
end module mpas_core
Modified: branches/omp_blocks/openmp_test/src/core_sw_elements/mpas_sw_time_integration.F
===================================================================
--- branches/omp_blocks/openmp_test/src/core_sw_elements/mpas_sw_time_integration.F 2012-08-23 21:56:35 UTC (rev 2122)
+++ branches/omp_blocks/openmp_test/src/core_sw_elements/mpas_sw_time_integration.F 2012-08-24 21:18:43 UTC (rev 2123)
@@ -36,17 +36,12 @@
stop
end if
- !$omp single private(block)
block => domain % blocklist
do while (associated(block))
block_d = block
- !$omp task firstprivate(block_d)
block_d % state % time_levs(2) % state % xtime % scalar = timeStamp
- !$omp end task
block => block % next
end do
- !$omp taskwait
- !$omp end single
end subroutine sw_timestep
@@ -76,10 +71,8 @@
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
- !$omp single
call mpas_timer_start('computations')
call mpas_setup_provis_states(domain % blocklist)
- !$omp end single
!
! Initialize time_levs(2) with state at current time
@@ -90,43 +83,31 @@
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(iCell, k)
do iCell=1,block % mesh % nCells ! couple tracers to h
do k=1,block % mesh % nVertLevels
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
- !$omp sections
- !$omp section
rk_weights(1) = dt/6.
rk_weights(2) = dt/3.
rk_weights(3) = dt/3.
rk_weights(4) = dt/6.
- !$omp section
rk_substep_weights(1) = dt/2.
rk_substep_weights(2) = dt/2.
rk_substep_weights(3) = dt
rk_substep_weights(4) = 0.
- !$omp end sections
- !$omp single
call mpas_timer_stop('computations')
- !$omp end single
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN RK loop
@@ -135,7 +116,6 @@
! --- update halos for diagnostic variables
- !$omp single
call mpas_timer_start('communications')
call mpas_dmpar_exch_halo_field(domain % blocklist % provis % pv_edge)
@@ -145,35 +125,25 @@
end if
call mpas_timer_stop('communications')
call mpas_timer_start('computations')
- !$omp end single
! --- compute tendencies
block => domain % blocklist
do while (associated(block))
- !$omp single
call mpas_timer_start('sw_compute_tend')
- !$omp end single
call sw_compute_tend(block % tend, block % provis, block % mesh)
- !$omp single
call mpas_timer_stop('sw_compute_tend')
call mpas_timer_start('sw_compute_scalar_tend')
- !$omp end single
call sw_compute_scalar_tend(block % tend, block % provis, block % mesh)
- !$omp single
call mpas_timer_stop('sw_compute_scalar_tend')
call mpas_timer_start('sw_enforce_boundary')
- !$omp end single
call sw_enforce_boundary_edge(block % tend, block % mesh)
- !$omp single
call mpas_timer_stop('sw_enforce_boundary')
- !$omp end single
block => block % next
end do
! --- update halos for prognostic variables
- !$omp single
call mpas_timer_stop('computations')
call mpas_timer_start('communications')
call mpas_dmpar_exch_halo_field(domain % blocklist % tend % u)
@@ -181,21 +151,17 @@
call mpas_dmpar_exch_halo_field(domain % blocklist % tend % tracers)
call mpas_timer_stop('communications')
call mpas_timer_start('computations')
- !$omp end single
! --- compute next substep state
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(iCell, k)
do iCell=1,block % mesh % nCells
do k=1,block % mesh % nVertLevels
block % provis % tracers % array(:,k,iCell) = ( block % state % time_levs(1) % state % h % array(k,iCell) * &
@@ -204,19 +170,12 @@
) / block % provis % h % array(k,iCell)
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
- !$omp single
call mpas_timer_start('sw_compute_solve_diagnostics')
- !$omp end single
call sw_compute_solve_diagnostics(dt, block % provis, block % mesh)
- !$omp single
call mpas_timer_stop('sw_compute_solve_diagnostics')
- !$omp end single
block => block % next
end do
end if
@@ -225,14 +184,11 @@
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(iCell, k)
do iCell=1,block % mesh % nCells
do k=1,block % mesh % nVertLevels
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
@@ -240,13 +196,10 @@
+ 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('computations')
- !$omp end single
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! END RK loop
@@ -256,12 +209,9 @@
!
! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
!
- !$omp single
call mpas_timer_start('computations')
- !$omp end single
block => domain % blocklist
do while (associated(block))
- !$omp do private(iCell, k)
do iCell=1,block % mesh % nCells
do k=1,block % mesh % nVertLevels
block % state % time_levs(2) % state % tracers % array(:,k,iCell) = &
@@ -269,17 +219,13 @@
/ block % state % time_levs(2) % state % h % array(k,iCell)
end do
end do
- !$omp end do
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
call sw_compute_solve_diagnostics(dt, block % state % time_levs(2) % state, block % mesh)
- !$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, &
@@ -287,15 +233,12 @@
block % state % time_levs(2) % state % uReconstructZonal % array, &
block % state % time_levs(2) % state % uReconstructMeridional % array &
)
- !$omp end single
block => block % next
end do
- !$omp single
call mpas_deallocate_provis_states(domain % blocklist)
call mpas_timer_stop('computations')
- !$omp end single
end subroutine sw_rk4
@@ -326,11 +269,12 @@
circulation, vorticity, ke, pv_edge, divergence, h_vertex
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
+ integer, dimension(:,:), pointer :: edgeSignOnCell, edgeSignOnVertex
real (kind=RKIND) :: r, u_diffusion
- real (kind=RKIND), save, allocatable, dimension(:,:) :: delsq_divergence
- real (kind=RKIND), save, allocatable, dimension(:,:) :: delsq_u
- real (kind=RKIND), save, allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+ real (kind=RKIND), pointer, dimension(:,:) :: delsq_divergence
+ real (kind=RKIND), pointer, dimension(:,:) :: delsq_u
+ real (kind=RKIND), pointer, dimension(:,:) :: delsq_circulation, delsq_vorticity
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
@@ -346,6 +290,10 @@
ke => s % ke % array
pv_edge => s % pv_edge % array
vh => s % vh % array
+ delsq_divergence => s % delsq_divergence % array
+ delsq_u => s % delsq_u % array
+ delsq_circulation => s % delsq_circulation % array
+ delsq_vorticity => s % delsq_vorticity % array
weightsOnEdge => grid % weightsOnEdge % array
kiteAreasOnVertex => grid % kiteAreasOnVertex % array
@@ -364,6 +312,8 @@
h_s => grid % h_s % array
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+ edgeSignOnVertex => grid % edgeSignOnVertex % array
tend_h => tend % h % array
tend_u => tend % u % array
@@ -379,6 +329,8 @@
meshScalingDel2 => grid % meshScalingDel2 % array
meshScalingDel4 => grid % meshScalingDel4 % array
+ !$omp parallel
+
!
! Compute height tendency for each cell
!
@@ -386,51 +338,24 @@
tend_h(:,:) = 0.0
!$omp end workshare
- !DWJ 08/22/12 Old tend_h loop
- if(config_old_loops) then
- !$omp do private(iEdge, cell1, cell2, flux)
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- !$omp critical (crit_cell1_tend_h)
- tend_h(k,cell1) = tend_h(k,cell1) - flux
- !$omp end critical (crit_cell1_tend_h)
-
- !$omp critical (crit_cell2_tend_h)
- tend_h(k,cell2) = tend_h(k,cell2) + flux
- !$omp end critical (crit_cell2_tend_h)
- end do
- end do
- !$omp end do
- else
!DWJ 08/22/12 New tend_h loop
!$omp do private(iCell, i, iEdge, k, flux)
do iCell=1,grid % nCellsSolve
do i = 1, grid % nEdgesOnCell % array(iCell)
iEdge = grid % edgesOnCell % array(i, iCell)
do k = 1, nVertLevels
- if(iCell == cellsOnEdge(1, iEdge)) then
- flux = - u(k, iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- else
- flux = u(k, iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- end if
- tend_h(k, iCell) = tend_h(k, iCell) + flux
+! if(iCell == cellsOnEdge(1, iEdge)) then
+! flux = - u(k, iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+! else
+! flux = u(k, iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+! end if
+ flux = edgeSignOnCell(i, iCell) * u(k, iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
+ tend_h(k, iCell) = tend_h(k, iCell) + flux / areaCell(iCell)
end do
end do
end do
!$omp end do
- end if
- !$omp do private(iCell, k)
- do iCell=1,grid % nCellsSolve
- do k=1,nVertLevels
- tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
- end do
- end do
- !$omp end do
-
#ifdef LANL_FORMULATION
!
! Compute u (normal) velocity tendency for each edge (cell face)
@@ -439,7 +364,7 @@
tend_u(:,:) = 0.0
!$omp end workshare
- !$omp do private(iCell, cell1, cell2, vertex1, vertex2, k, q, j, eoe, workpv)
+ !$omp do private(iEdge, cell1, cell2, vertex1, vertex2, k, q, j, eoe, workpv)
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -474,7 +399,7 @@
tend_u(:,:) = 0.0
!$omp end workshare
- !$omp do private(iEdge, vertex1, vertex2, cell1, cell2, k, vorticity_abs, workpv)
+ !$omp do private(iEdge, vertex1, vertex2, cell1, cell2, k, vorticity_abs, workpv)
do iEdge=1,grid % nEdgesSolve
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
@@ -500,7 +425,7 @@
! Compute diffusion, computed as </font>
<font color="black">abla divergence - k \times </font>
<font color="gray">abla vorticity
! only valid for visc == constant
if (config_h_mom_eddy_visc2 > 0.0) then
- !$omp do private(iEdge, cell1, cell2, vertex1, vertex2, k, u_diffusion)
+ !$omp do private(iEdge, cell1, cell2, vertex1, vertex2, k, u_diffusion)
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -524,23 +449,12 @@
! strictly only valid for h_mom_eddy_visc4 == constant
!
if (config_h_mom_eddy_visc4 > 0.0) then
- !$omp sections
- !$omp section
- allocate(delsq_divergence(nVertLevels, nCells+1))
- !$omp section
- allocate(delsq_u(nVertLevels, nEdges+1))
- !$omp section
- allocate(delsq_circulation(nVertLevels, nVertices+1))
- !$omp section
- allocate(delsq_vorticity(nVertLevels, nVertices+1))
- !$omp end sections
-
!$omp workshare
delsq_u(:,:) = 0.0
!$omp end workshare
! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity
- !$omp do private(iEdge, cell1, cell2, vertex1, vertex2, k)
+ !$omp do private(iEdge, cell1, cell2, vertex1, vertex2, k)
do iEdge=1,grid % nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -561,34 +475,6 @@
delsq_circulation(:,:) = 0.0
!$omp end workshare
- if(config_old_loops) then
- !DWJ 08/22/12 Old delsq_circ loop
- !$omp do private(iEdge, vertex1, vertex2)
- do iEdge=1,nEdges
- vertex1 = verticesOnEdge(1,iEdge)
- vertex2 = verticesOnEdge(2,iEdge)
- do k=1,nVertLevels
- !$omp critical (crit_vert1_delsq_circ)
- delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
- - dcEdge(iEdge) * delsq_u(k,iEdge)
- !$omp end critical (crit_vert1_delsq_circ)
-
- !$omp critical (crit_vert2_delsq_cric)
- delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
- + dcEdge(iEdge) * delsq_u(k,iEdge)
- !$omp end critical (crit_vert2_delsq_cric)
- end do
- end do
- !$omp end do
- !$omp do private(iVertex, r, k)
- do iVertex=1,nVertices
- r = 1.0 / areaTriangle(iVertex)
- do k=1,nVertLevels
- delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
- end do
- end do
- !$omp end do
- else
!DWJ 08/22/12 New delsq_circ loop
!$omp do private(iVertex, i, iEdge, k)
do iVertex=1,nVertices
@@ -600,62 +486,33 @@
end do
end do
!$omp end do
- end if
! Divergence using </font>
<font color="red">abla^2 u
!$omp workshare
delsq_divergence(:,:) = 0.0
!$omp end workshare
- if(config_old_loops) then
- !DWJ 08/22/12 Old delsq_div loop
- !$omp do private(iEdge, cell1, cell2, k)
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,nVertLevels
- !$omp critical (crit_cell1_delsq_div)
- delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
- + delsq_u(k,iEdge)*dvEdge(iEdge)
- !$omp end critical (crit_cell1_delsq_div)
-
- !$omp critical (crit_cell2_delsq_div)
- delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
- - delsq_u(k,iEdge)*dvEdge(iEdge)
- !$omp end critical (crit_cell2_delsq_div)
- end do
- end do
- !$omp end do
-
- !$omp do private(iCell, r, k)
- do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
- do k = 1,nVertLevels
- delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
- end do
- end do
- !$omp end do
- else
!DWJ 08/22/12 New delsq_div loop
!$omp do private(iCell, i, iEdge, k)
do iCell = 1, nCells
do i = 1, nEdgesOnCell(iCell)
iEdge = edgesOnEdge(i, iCell)
do k = 1, nVertLevels
- if(iCell == cellsOnEdge(1, iEdge)) then
- delsq_divergence(k, iCell) = delsq_divergence(k, iCell) + delsq_u(k, iEdge)*dvEdge(iEdge) / areaCell(iCell)
- else
- delsq_divergence(k, iCell) = delsq_divergence(k, iCell) - delsq_u(k, iEdge)*dvEdge(iEdge) / areaCell(iCell)
- end if
+! if(iCell == cellsOnEdge(1, iEdge)) then
+! delsq_divergence(k, iCell) = delsq_divergence(k, iCell) + delsq_u(k, iEdge)*dvEdge(iEdge) / areaCell(iCell)
+! else
+! delsq_divergence(k, iCell) = delsq_divergence(k, iCell) - delsq_u(k, iEdge)*dvEdge(iEdge) / areaCell(iCell)
+! end if
+
+ delsq_divergence(k, iCell) = delsq_divergence(k, iCell) + edgeSignOnCell(i, iCell) * delsq_u(k, iEdge)*dvEdge(iEdge) / areaCell(iCell)
end do
end do
end do
!$omp end do
- end if
! 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(iEdge, cell1, cell2, vertex1, vertex2, k, u_diffusion)
+ !$omp do private(iEdge, cell1, cell2, vertex1, vertex2, k, u_diffusion)
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -675,23 +532,11 @@
end do
end do
!$omp end do
-
- !$omp sections
- !$omp section
- deallocate(delsq_divergence)
- !$omp section
- deallocate(delsq_u)
- !$omp section
- deallocate(delsq_circulation)
- !$omp section
- deallocate(delsq_vorticity)
- !$omp end sections
-
end if
! Compute u (velocity) tendency from wind stress (u_src)
if(config_wind_stress) then
- !$omp do private(iEdge)
+ !$omp do private(iEdge)
do iEdge=1,grid % nEdges
tend_u(1,iEdge) = tend_u(1,iEdge) &
+ u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
@@ -700,7 +545,7 @@
endif
if (config_bottom_drag) then
- !$omp do private(iEdge, ke_edge)
+ !$omp do private(iEdge, ke_edge)
do iEdge=1,grid % nEdges
! bottom drag is the same as POP:
! -c |u| u where c is unitless and 1.0e-3.
@@ -714,6 +559,7 @@
end do
!$omp end do
endif
+ !$omp end parallel
end subroutine sw_compute_tend
@@ -737,13 +583,14 @@
real (kind=RKIND) :: flux, tracer_edge, r
real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux
integer, dimension(:,:), pointer :: boundaryEdge
- real (kind=RKIND), dimension(:,:), save, allocatable :: boundaryMask
- real (kind=RKIND), dimension(:,:,:), save, allocatable:: delsq_tracer
+
+ real (kind=RKIND), dimension(:,:), pointer :: boundaryMask
+ real (kind=RKIND), dimension(:,:,:), pointer:: delsq_tracer
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
real (kind=RKIND), dimension(:,:,:), pointer :: tracers, tracer_tend
- integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell
+ integer, dimension(:,:), pointer :: cellsOnEdge, boundaryCell, edgeSignOnCell, edgeSignOnVertex
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
real (kind=RKIND), dimension(:,:), pointer :: u, h_edge
@@ -759,40 +606,22 @@
boundaryEdge=> grid % boundaryEdge % array
areaCell => grid % areaCell % array
tracer_tend => tend % tracers % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+ edgeSignOnVertex => grid % edgeSignOnVertex % array
+ boundaryMask => grid % boundaryMask % array
+ delsq_tracer => s % delsq_tracer % array
+
coef_3rd_order = 0.
if (config_tracer_adv_order == 3) coef_3rd_order = 1.0
if (config_tracer_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
- !$omp workshare
+ !$omp parallel
+
tracer_tend(:,:,:) = 0.0
- !$omp end workshare
if (config_tracer_adv_order == 2) then
- if(config_old_loops) then
- !$omp do private(iEdge, cell1, cell2, k, iTracer, tracer_edge, flux)
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
- do k=1,grid % nVertLevels
- do iTracer=1,grid % nTracers
- tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
- !$omp critical (crit_cell1_2nd_tracer_tend_hadv)
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- !$omp end critical (crit_cell1_2nd_tracer_tend_hadv)
-
- !$omp critical (crit_cell2_2nd_tracer_tend_hadv)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
- !$omp end critical (crit_cell2_2nd_tracer_tend_hadv)
- end do
- end do
- end if
- end do
- !$omp end do
- else
!$omp do private(iCell, i, iEdge, cell1, cell2, k, iTracer, tracer_edge, flux)
do iCell = 1, grid % nCellsSolve
do i = 1, grid % nEdgesOnCell % array(iCell)
@@ -802,12 +631,14 @@
do k = 1, grid % nVertLevels
do iTracer = 1, grid % nTracers
tracer_edge= 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
- if(iCell == cell1) then
- flux = - u(k, iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
- else
- flux = u(k, iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
- end if
+! if(iCell == cell1) then
+! flux = - u(k, iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
+! else
+! flux = u(k, iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
+! end if
+ flux = edgeSignOnCell(i, iCell) * u(k, iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
+
tracer_tend(iTracer, k, iCell) = tracer_tend(iTracer,k,iCell) + flux/areaCell(iCell)
end do
@@ -815,74 +646,9 @@
end do
end do
!$omp end do
- end if
else if (config_tracer_adv_order == 3) then
- if(config_old_loops) then
- !$omp do private(iEdge, cell1, cell2, d2fdx2_cell1, d2fdx2_cell2, iTracer, i, flux)
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if a cell not on the most outside ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- do iTracer=1,grid % nTracers
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- !-- else u <= 0:
- else
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
- end if
-
- !-- update tendency
- !$omp critical (crit_cell1_3rd_tracer_tend_hadv)
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- !$omp end critical (crit_cell1_3rd_tracer_tend_hadv)
-
- !$omp critical (crit_cell2_3rd_tracer_tend_hadv)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
- !$omp end critical (crit_cell2_3rd_tracer_tend_hadv)
- enddo
- end do
- end if
- end do
- !$omp end do
- else
!$omp do private(iCell, j, iEdge, cell1, cell2, k, d2fdx2_cell1, d2fdx2_cell2, iTracer, i, flux)
do iCell = 1, grid % nCellsSolve
do j = 1, grid % nEdgesOnCell % array(iCell)
@@ -930,75 +696,22 @@
+(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
end if
- if(cell1 == iCell) then
- tracer_tend(iTracer,k,iCell) = tracer_tend(iTracer,k,iCell) - flux / areaCell(iCell)
- else
- tracer_tend(iTracer,k,iCell) = tracer_tend(iTracer,k,iCell) + flux / areaCell(iCell)
- end if
+! if(cell1 == iCell) then
+! tracer_tend(iTracer,k,iCell) = tracer_tend(iTracer,k,iCell) - flux / areaCell(iCell)
+! else
+! tracer_tend(iTracer,k,iCell) = tracer_tend(iTracer,k,iCell) + flux / areaCell(iCell)
+! end if
+
+ tracer_tend(iTracer,k,iCell) = tracer_tend(iTracer,k,iCell) + edgeSignOncell(i, iCell) * flux / areaCell(iCell)
enddo
end do
end do
end do
!$omp end do
- end if
else if (config_tracer_adv_order == 4) then
- if(config_old_loops) then
- !$omp do private(iEdge, cell1, cell2, k, d2fdx2_cell1, d2fdx2_cell2, iTracer, i, flux)
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- !-- if an edge is not on the outer-most ring of the halo
- if (cell1 <= grid%nCells .and. cell2 <= grid%nCells) then
-
- do k=1,grid % nVertLevels
-
- d2fdx2_cell1 = 0.0
- d2fdx2_cell2 = 0.0
-
- do iTracer=1,grid % nTracers
-
- !-- if not a boundary cell
- if(boundaryCell(k,cell1).eq.0.and.boundaryCell(k,cell2).eq.0) then
-
- d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracers(iTracer,k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracers(iTracer,k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell1))
- end do
-
- !-- all edges of cell 2
- do i=1, grid % nEdgesOnCell % array (cell2)
- d2fdx2_cell2 = d2fdx2_cell2 + &
- deriv_two(i+1,2,iEdge) * tracers(iTracer,k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- flux = dvEdge(iEdge) * u(k,iEdge) * h_edge(k,iEdge) * ( &
- 0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
-
- !-- update tendency
- !$omp critical (crit_cell1_4th_tracer_tend_hadv)
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- !$omp end critical (crit_cell1_4th_tracer_tend_hadv)
-
- !$omp critical (crit_cell2_4th_tracer_tend_hadv)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
- !$omp end critical (crit_cell2_4th_tracer_tend_hadv)
- enddo
- end do
- end if
- end do
- !$omp end do
- else
- !$omp do private(iCell, j, iEdge, cell1, cell2, k, d2fdx2_cell1, d2fdx2_cell2, iTracer, i)
+ !$omp do private(iCell, j, iEdge, cell1, cell2, k, d2fdx2_cell1, d2fdx2_cell2, iTracer, i, flux)
do iCell = 1, grid % nCellsSolve
do j = 1, grid % nEdgesOnCell % array(iCell)
iEdge = grid % edgesOnCell % array(j, iCell)
@@ -1035,18 +748,19 @@
0.5*(tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2)) &
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
- if(cell1 == iCell) then
- tracer_tend(iTracer, k, iCell) = tracer_tend(iTracer, k, iCell) - flux/areaCell(iCell)
- else
- tracer_tend(iTracer, k, iCell) = tracer_tend(iTracer, k, iCell) + flux/areaCell(iCell)
- end if
+! if(cell1 == iCell) then
+! tracer_tend(iTracer, k, iCell) = tracer_tend(iTracer, k, iCell) - flux/areaCell(iCell)
+! else
+! tracer_tend(iTracer, k, iCell) = tracer_tend(iTracer, k, iCell) + flux/areaCell(iCell)
+! end if
+
+ tracer_tend(iTracer, k, iCell) = tracer_tend(iTracer, k, iCell) + edgeSignOnCell(i, iCell) * flux/areaCell(iCell)
enddo
end do
end do
end do
!$omp end do
- end if
endif ! if (config_tracer_adv_order == 2 )
@@ -1058,74 +772,37 @@
!
! compute a boundary mask to enforce insulating boundary conditions in the horizontal
!
- !$omp single
- allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
- !$omp end single
-
!$omp workshare
boundaryMask = 1.0
where(boundaryEdge.eq.1) boundaryMask=0.0
!$omp end workshare
- if(config_old_loops) then
- !$omp do private(iEdge, cell1, cell2, invAreaCell1, invAreaCell2, k, iTracer, tracer_turb_flux, flux)
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- invAreaCell1 = 1.0/areaCell(cell1)
- invAreaCell2 = 1.0/areaCell(cell2)
+ !$omp do private(iCell, invAreaCell1, i, iEdge, cell1, cell2, k, iTracer, tracer_turb_flux, flux)
+ do iCell = 1, grid % nCellsSolve
+ invAreaCell1 = 1.0 / areaCell(iCell)
+ do i = 1, grid % nEdgesOnCell % array(iCell)
+ iEdge = grid % edgesOnCell % array(i, iCell)
+ cell1 = cellsOnEdge(1, iEdge)
+ cell2 = cellsOnEdge(2, iEdge)
- do k=1,grid % nVertLevels
- do iTracer=1, grid % nTracers
- ! \kappa_2 </font>
<font color="red">abla \phi on edge
- tracer_turb_flux = config_h_tracer_eddy_diff2 &
- *( tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge)
+ do k = 1, grid %nVertLevels
+ do iTracer = 1, grid % nTracers
+ tracer_turb_flux = config_h_tracer_eddy_diff2 * (tracers(iTracer,k, cell2) - tracers(iTracer, k, cell1)) /dcEdge(iEdge)
- ! div(h \kappa_2 </font>
<font color="gray">abla \phi) at cell center
- flux = dvEdge(iEdge) * h_edge(k,iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
- !$omp critical (crit_cell1_tracer_tend_del2)
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) + flux * invAreaCell1
- !$omp end critical (crit_cell1_tracer_tend_del2)
+ flux = dvEdge(iEdge) * h_edge(k, iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
- !$omp critical (crit_cell2_tracer_tend_del2)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) - flux * invAreaCell2
- !$omp end critical (crit_cell2_tracer_tend_del2)
- end do
- end do
+! if(iCell == cell1) then
+! tracer_tend(iTracer, k, iCell) = tracer_tend(iTracer, k, iCell) + flux * invAreaCell1
+! else
+! tracer_tend(iTracer, k, iCell) = tracer_tend(iTracer, k, iCell) - flux * invAreaCell1
+! end if
- end do
- !$omp end do
- else
- !$omp do private(iCell, invAreaCell1, i, iEdge, cell1, cell2, k, iTracer, tracer_turb_flux, flux)
- do iCell = 1, grid % nCellsSolve
- invAreaCell1 = 1.0 / areaCell(iCell)
- do i = 1, grid % nEdgesOnCell % array(iCell)
- iEdge = grid % edgesOnCell % array(i, iCell)
- cell1 = cellsOnEdge(1, iEdge)
- cell2 = cellsOnEdge(2, iEdge)
-
- do k = 1, grid %nVertLevels
- do iTracer = 1, grid % nTracers
- tracer_turb_flux = config_h_tracer_eddy_diff2 * (tracers(iTracer,k, cell2) - tracers(iTracer, k, cell1)) /dcEdge(iEdge)
-
- flux = dvEdge(iEdge) * h_edge(k, iEdge) * tracer_turb_flux * boundaryMask(k, iEdge)
-
- if(iCell == cell1) then
- tracer_tend(iTracer, k, iCell) = tracer_tend(iTracer, k, iCell) + flux * invAreaCell1
- else
- tracer_tend(iTracer, k, iCell) = tracer_tend(iTracer, k, iCell) - flux * invAreaCell1
- end if
- end do
+ tracer_tend(iTracer, k, iCell) = tracer_tend(iTracer, k, iCell) - edgeSignOnCell(i,iCell) * flux * invAreaCell1
end do
end do
end do
- !$omp end do
- end if
-
- !$omp single
- deallocate(boundaryMask)
- !$omp end single
-
+ end do
+ !$omp end do
end if
!
@@ -1137,132 +814,78 @@
!
! compute a boundary mask to enforce insulating boundary conditions in the horizontal
!
- !$omp sections
- !$omp section
- allocate(boundaryMask(grid % nVertLevels, grid % nEdges+1))
- !$omp section
- allocate(delsq_tracer(grid % nTracers, grid % nVertLevels, grid % nCells+1))
- !$omp end sections
-
!$omp workshare
boundaryMask = 1.0
where(boundaryEdge.eq.1) boundaryMask=0.0
delsq_tracer(:,:,:) = 0.
!$omp end workshare
- if(config_old_loops) then
- ! first del2: div(h </font>
<font color="red">abla \phi) at cell center
- !$omp do private(iEdge, cell1, cell2, k, iTracer)
- do iEdge=1,grid % nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ !DWJ 08/23/12 New Del4 loop 1
+ !$omp do private(iCell, i, iEdge, cell1, cell2, k, iTracer)
+ do iCell = 1, grid % nCellsSolve
+ do i = 1, grid % nEdgesOnCell % array(iCell)
+ iEdge = grid % edgesOnCell % array(i, iCell)
+ cell1 = cellsOnEdge(1, iEdge)
+ cell2 = cellsOnEdge(2, iEdge)
- do k=1,grid % nVertLevels
- do iTracer=1, grid % nTracers
- delsq_tracer(iTracer,k,cell1) = delsq_tracer(iTracer,k,cell1) &
- + dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
- delsq_tracer(iTracer,k,cell2) = delsq_tracer(iTracer,k,cell2) &
- - dvEdge(iEdge) * h_edge(k,iEdge) * (tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) / dcEdge(iEdge) * boundaryMask(k,iEdge)
- end do
- end do
+ do k = 1, grid % nVertLevels
+ do iTracer = 1, grid % nTracers
+! if(iCell == cell1) then
+! delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer, k, iCell) + dvEdge(iEdge) * h_edge(k,iEdge) &
+! * (tracers(iTracer, k, cell2) - tracers(iTracer, k, cell1))/dcEdge(iEdge) * boundaryMask(k,iEdge)
+! else
+! delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer, k, iCell) - dvEdge(iEdge) * h_edge(k,iEdge) &
+! * (tracers(iTracer, k, cell2) - tracers(iTracer, k, cell1))/dcEdge(iEdge) * boundaryMask(k,iEdge)
+! end if
- end do
- !$omp end do
- else
- !$omp do private(iCell, i, iEdge, cell1, cell2, k, iTracer)
- do iCell = 1, grid % nCellsSolve
- do i = 1, grid % nEdgesOnCell % array(iCell)
- iEdge = grid % edgesOnCell % array(i, iCell)
- cell1 = cellsOnEdge(1, iEdge)
- cell2 = cellsOnEdge(2, iEdge)
-
- do k = 1, grid % nVertLevels
- do iTracer = 1, grid % nTracers
- if(iCell == cell1) then
- delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer, k, iCell) + dvEdge(iEdge) * h_edge(k,iEdge) &
- * (tracers(iTracer, k, cell2) - tracers(iTracer, k, cell1))/dcEdge(iEdge) * boundaryMask(k,iEdge)
- else
- delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer, k, iCell) - dvEdge(iEdge) * h_edge(k,iEdge) &
- * (tracers(iTracer, k, cell2) - tracers(iTracer, k, cell1))/dcEdge(iEdge) * boundaryMask(k,iEdge)
- end if
- end do
+ delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer, k, iCell) - edgeSignOnCell(i, iCell) * dvEdge(iEdge) * h_edge(k,iEdge) &
+ * (tracers(iTracer, k, cell2) - tracers(iTracer, k, cell1))/dcEdge(iEdge) * boundaryMask(k,iEdge)
end do
end do
end do
- !$omp end do
- end if
+ end do
+ !$omp end do
- !$omp do private(iCell, r, k, iTracer)
+ !$omp do private(iCell, r, k, iTracer)
do iCell = 1, grid % nCells
r = 1.0 / grid % areaCell % array(iCell)
do k=1,grid % nVertLevels
- do iTracer=1,grid % nTracers
- delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
+ do iTracer=1,grid % nTracers
+ delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
+ end do
end do
- end do
end do
!$omp end do
- if(config_old_loops) then
- ! second del2: div(h </font>
<font color="gray">abla [delsq_tracer]) at cell center
- !$omp do private(iEdge, cell1, cell2, invAreaCell1, invAreaCell2, k, iTracer, tracer_turb_flux, flux)
- do iEdge=1,grid % nEdges
- cell1 = grid % cellsOnEdge % array(1,iEdge)
- cell2 = grid % cellsOnEdge % array(2,iEdge)
- invAreaCell1 = 1.0 / grid % areaCell % array(cell1)
- invAreaCell2 = 1.0 / grid % areaCell % array(cell2)
+ !$omp do private(iCell, invAreaCell1, i, iEdge, cell1, cell2, k, iTracer, tracer_turb_flux, flux)
+ do iCell = 1, grid % nCellsSolve
+ invAreaCell1 = 1.0 / grid % areaCell % array(iCell)
+ do i = 1, grid % nEdgesOnCell % array(iCell)
+ iEdge = grid % edgesOnCell % array(i, iCell)
+ cell1 = grid % cellsOnEdge % array(1, iEdge)
+ cell2 = grid % cellsOnEdge % array(2, iEdge)
+ do k = 1, grid % nVertLevels
+ do iTracer = 1, grid % nTracers
+ tracer_turb_flux = config_h_tracer_eddy_diff4 * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) /dcEdge(iEdge)
+ flux = dvEdge(iEdge) * tracer_turb_flux
- do k=1,grid % nVertLevels
- do iTracer=1,grid % nTracers
- tracer_turb_flux = config_h_tracer_eddy_diff4 * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) / dcEdge(iEdge)
- flux = dvEdge(iEdge) * tracer_turb_flux
- !$omp critical (crit_cell1_tracer_tend_del4)
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux * invAreaCell1 * boundaryMask(k,iEdge)
- !$omp end critical (crit_cell1_tracer_tend_del4)
+! if(iCell == cell1) then
+! tracer_tend(iTracer,k,iCell) = tracer_tend(iTracer,k,iCell) - flux * invAreaCell1 * boundaryMask(k,iEdge)
+! else
+! tracer_tend(iTracer,k,iCell) = tracer_tend(iTracer,k,iCell) + flux * invAreaCell1 * boundaryMask(k,iEdge)
+! end if
- !$omp critical (crit_cell2_tracer_tend_del4)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux * invAreaCell2 * boundaryMask(k,iEdge)
- !$omp end critical (crit_cell2_tracer_tend_del4)
- end do
- enddo
-
- end do
- !$omp end do
- else
- !$omp do private(iCell, invAreaCell1, i, iEdge, cell1, cell2, k, iTracer, tracer_turb_flux, flux)
- do iCell = 1, grid % nCellsSolve
- invAreaCell1 = 1.0 / grid % areaCell % array(iCell)
- do i = 1, grid % nEdgesOnCell % array(iCell)
- iEdge = grid % edgesOnCell % array(i, iCell)
- cell1 = grid % cellsOnEdge % array(1, iEdge)
- cell2 = grid % cellsOnEdge % array(2, iEdge)
- do k = 1, grid % nVertLevels
- do iTracer = 1, grid % nTracers
- tracer_turb_flux = config_h_tracer_eddy_diff4 * (delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)) /dcEdge(iEdge)
- flux = dvEdge(iEdge) * tracer_turb_flux
-
- if(iCell == cell1) then
- tracer_tend(iTracer,k,iCell) = tracer_tend(iTracer,k,iCell) - flux * invAreaCell1 * boundaryMask(k,iEdge)
- else
- tracer_tend(iTracer,k,iCell) = tracer_tend(iTracer,k,iCell) + flux * invAreaCell1 * boundaryMask(k,iEdge)
- end if
- end do
+ tracer_tend(iTracer,k,iCell) = tracer_tend(iTracer,k,iCell) + edgeSignOnCell(i, iCell) * flux * invAreaCell1 * boundaryMask(k,iEdge)
end do
end do
end do
- !$omp end do
- end if
+ end do
+ !$omp end do
+ end if
- !$omp sections
- !$omp section
- deallocate(delsq_tracer)
- !$omp section
- deallocate(boundaryMask)
- !$omp end sections
+ !$omp end parallel
- end if
-
end subroutine sw_compute_scalar_tend
@@ -1291,6 +914,7 @@
circulation, vorticity, ke, pv_edge, pv_vertex, pv_cell, gradPVn, gradPVt, divergence, &
h_vertex, vorticity_cell
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex, boundaryEdge, boundaryCell
+ integer, dimension(:,:), pointer :: edgeSignOnCell, edgeSignOnVertex
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, h1, h2
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
@@ -1334,6 +958,8 @@
fVertex => grid % fVertex % array
fEdge => grid % fEdge % array
deriv_two => grid % deriv_two % array
+ edgeSignOnCell => grid % edgeSignOnCell % array
+ edgeSignOnVertex => grid % edgeSignOnVertex % array
nCells = grid % nCells
nEdges = grid % nEdges
@@ -1343,6 +969,8 @@
boundaryEdge => grid % boundaryEdge % array
boundaryCell => grid % boundaryCell % array
+ !$omp parallel
+
!
! Find those cells that have an edge on the boundary
!
@@ -1350,36 +978,17 @@
boundaryCell(:,:) = 0
!$omp end workshare
- if(config_old_loops) then
- !$omp do private(iEdge, k, cell1, cell2)
- do iEdge=1,nEdges
- do k=1,nVertLevels
- if(boundaryEdge(k,iEdge).eq.1) then
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- !$omp critical (crit_cell1_boundary_cell)
- boundaryCell(k,cell1) = 1
- !$omp end critical (crit_cell1_boundary_cell)
- !$omp critical (crit_cell2_boundary_cell)
- boundaryCell(k,cell2) = 1
- !$omp end critical (crit_cell2_boundary_cell)
- endif
- enddo
- enddo
- !$omp end do
- else
- !$omp do private(iCell, i, iEdge, k)
- do iCell = 1, nCells
- do i = 1, grid % nEdgesOnCell % array(iCell)
- iEdge = grid % edgesOnCell % array(i, iCell)
- do k = 1, nVertLevels
- boundaryCell(k, iCell) = max(boundaryCell(k,iCell), boundaryEdge(k, iEdge))
- end do
+ !$omp do private(iCell, i, iEdge, k)
+ do iCell = 1, nCells
+ do i = 1, grid % nEdgesOnCell % array(iCell)
+ iEdge = grid % edgesOnCell % array(i, iCell)
+ do k = 1, nVertLevels
+ boundaryCell(k, iCell) = max(boundaryCell(k,iCell), boundaryEdge(k, iEdge))
end do
end do
- !$omp end do
- end if
+ end do
+ !$omp end do
!
! Compute height on cell edges at velocity locations
@@ -1392,7 +1001,7 @@
if (config_thickness_adv_order == 2) then
- !$omp do private(iEdge, cell1, cell2, k)
+ !$omp do private(iEdge, cell1, cell2, k)
do iEdge=1,grid % nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -1406,7 +1015,7 @@
else if (config_thickness_adv_order == 3) then
- !$omp do private(iEdge, cell1, cell2, k, d2fdx2_cell1, d2fdx2_cell2, i)
+ !$omp do private(iEdge, cell1, cell2, d2fdx2_cell1, d2fdx2_cell2, i)
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -1460,7 +1069,7 @@
else if (config_thickness_adv_order == 4) then
- !$omp do private(iEdge, cell1, cell2, k, d2fdx2_cell1, d2fdx2_cell2, i)
+ !$omp do private(iEdge, cell1, cell2, k, d2fdx2_cell1, d2fdx2_cell2, i)
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
@@ -1517,41 +1126,27 @@
circulation(:,:) = 0.0
!$omp end workshare
- if(config_old_loops) then
- !$omp do private(iEdge, k)
- do iEdge=1,nEdges
- do k=1,nVertLevels
- !$omp critical (crit_vert1_circ)
- circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
- !$omp end critical (crit_vert1_circ)
+ !$omp do private(iVertex, i, iEdge, vertex1, vertex2, k)
+ do iVertex = 1, grid % nVertices
+ do i = 1, grid % vertexDegree
+ iEdge = grid % edgesOnVertex % array(i, iVertex)
+ vertex1 = grid % verticesOnEdge % array(1, iEdge)
+ vertex2 = grid % verticesOnEdge % array(2, iEdge)
- !$omp critical (crit_vert2_circ)
- circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
- !$omp end critical (crit_vert2_circ)
- end do
- end do
- !$omp end do
- else
- !$omp do private(iVertex, i, iEdge, vertex1, vertex2, k)
- do iVertex = 1, grid % nVertices
- do i = 1, grid % vertexDegree
- iEdge = grid % edgesOnVertex % array(i, iVertex)
- vertex1 = grid % verticesOnEdge % array(1, iEdge)
- vertex2 = grid % verticesOnEdge % array(2, iEdge)
+ do k = 1, grid % nVertLevels
+! if(vertex1 == iVertex) then
+! circulation(k, iVertex) = circulation(k, iVertex) - dcEdge(iEdge) * u(k, iEdge)
+! else
+! circulation(k, iVertex) = circulation(k, iVertex) + dcEdge(iEdge) * u(k, iEdge)
+! end if
- do k = 1, grid % nVertLevels
- if(vertex1 == iVertex) then
- circulation(k, iVertex) = circulation(k, iVertex) - dcEdge(iEdge) * u(k, iEdge)
- else
- circulation(k, iVertex) = circulation(k, iVertex) + dcEdge(iEdge) * u(k, iEdge)
- end if
- end do
+ circulation(k, iVertex) = circulation(k, iVertex) + edgeSignOnVertex(i, iVertex) * dcEdge(iEdge) * u(k, iEdge)
end do
end do
- !$omp end do
- end if
+ end do
+ !$omp end do
- !$omp do private(iVertex, k)
+ !$omp do private(iVertex, k)
do iVertex=1,nVertices
do k=1,nVertLevels
vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
@@ -1567,48 +1162,27 @@
divergence(:,:) = 0.0
!$omp end workshare
- if(config_old_loops) then
- !$omp do private(iEdge, cell1, cell2, k)
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,nVertLevels
- !$omp critical (crit_cell1_div)
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
- !$omp end critical (crit_cell1_div)
- enddo
- endif
- if(cell2 <= nCells) then
- do k=1,nVertLevels
- !$omp critical (crit_cell2_div)
- divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
- !$omp end critical (crit_cell2_div)
- enddo
- end if
- end do
- !$omp end do
- else
- !$omp do private(iCell, i, iEdge, cell1, cell2)
- do iCell = 1, nCells
- do i = 1, grid % nEdgesOnCell % array(iCell)
- iEdge = grid % edgesOnCell % array(i, iCell)
- cell1 = cellsOnEdge(1, iEdge)
- cell2 = cellsOnEdge(2, iEdge)
+ !$omp do private(iCell, i, iEdge, cell1, cell2, k)
+ do iCell = 1, nCells
+ do i = 1, grid % nEdgesOnCell % array(iCell)
+ iEdge = grid % edgesOnCell % array(i, iCell)
+ cell1 = cellsOnEdge(1, iEdge)
+ cell2 = cellsOnEdge(2, iEdge)
- do k = 1, nVertLevels
- if(cell1 == iCell) then
- divergence(k, iCell) = divergence(k, iCell) + u(k, iEdge) * dvEdge(iEdge)
- else
- divergence(k, iCell) = divergence(k, iCell) - u(k, iEdge) * dvEdge(iEdge)
- end if
- end do
+ do k = 1, nVertLevels
+! if(cell1 == iCell) then
+! divergence(k, iCell) = divergence(k, iCell) + u(k, iEdge) * dvEdge(iEdge)
+! else
+! divergence(k, iCell) = divergence(k, iCell) - u(k, iEdge) * dvEdge(iEdge)
+! end if
+
+ divergence(k, iCell) = divergence(k, iCell) - edgeSignOnCell(i, iCell) * u(k, iEdge) * dvEdge(iEdge)
end do
end do
- !$omp end do
- end if
+ end do
+ !$omp end do
- !$omp do private(iCell, r, k)
+ !$omp do private(iCell, r, k)
do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
do k = 1,nVertLevels
@@ -1624,7 +1198,7 @@
ke(:,:) = 0.0
!$omp end workshare
- !$omp do private(iCell, i, iEdge, k)
+ !$omp do private(iCell, i, iEdge, k)
do iCell=1,nCells
do i=1,nEdgesOnCell(iCell)
iEdge = edgesOnCell(i,iCell)
@@ -1645,7 +1219,7 @@
v(:,:) = 0.0
!$omp end workshare
- !$omp do private(iEdge, i, eoe, k)
+ !$omp do private(iEdge, i, eoe, k)
do iEdge = 1,nEdges
do i=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(i,iEdge)
@@ -1664,7 +1238,7 @@
vh(:,:) = 0.0
!$omp end workshare
- !$omp do private(iEdge, j, eoe, k)
+ !$omp do private(iEdge, j, eoe, k)
do iEdge=1,grid % nEdgesSolve
do j=1,nEdgesOnEdge(iEdge)
eoe = edgesOnEdge(j,iEdge)
@@ -1676,12 +1250,11 @@
!$omp end do
#endif
-
!
! Compute height at vertices, pv at vertices, and average pv to edge locations
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
- !$omp do private(iVertex, k, i)
+ !$omp do private(iVertex, k, i)
do iVertex = 1,nVertices
do k=1,nVertLevels
h_vertex(k,iVertex) = 0.0
@@ -1699,7 +1272,7 @@
! Compute gradient of PV in the tangent direction
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
!
- !$omp do private(iEdge, k)
+ !$omp do private(iEdge, k)
do iEdge = 1,nEdges
do k = 1,nVertLevels
gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
@@ -1716,22 +1289,20 @@
pv_edge(:,:) = 0.0
!$omp end workshare
- !$omp do private(iVertex, i, iEdge, k)
- do iVertex = 1,nVertices
- do i=1,grid % vertexDegree
- iEdge = edgesOnVertex(i,iVertex)
- do k=1,nVertLevels
- pv_edge(k,iEdge) = pv_edge(k,iEdge) + 0.5 * pv_vertex(k,iVertex)
- end do
+ !$omp do private(iEdge, vertex1, vertex2, k)
+ do iEdge = 1, nEdges
+ vertex1 = verticesOnEdge(1, iEdge)
+ vertex2 = verticesOnEdge(2, iEdge)
+ do k = 1, nVertLevels
+ pv_edge(k, iEdge) = 0.5 * (pv_vertex(k, vertex1) + pv_vertex(k, vertex2))
end do
end do
!$omp end do
-
!
! Modify PV edge with upstream bias.
!
- !$omp do private(iEdge, k)
+ !$omp do private(iEdge, k)
do iEdge = 1,nEdges
do k = 1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * v(k,iEdge) * dt * gradPVt(k,iEdge)
@@ -1739,7 +1310,6 @@
enddo
!$omp end do
-
!
! Compute pv at cell centers
! ( this computes pv_cell for all real cells and distance-1 ghost cells )
@@ -1749,18 +1319,20 @@
vorticity_cell(:,:) = 0.0
!$omp end workshare
- !$omp do private(iVertex, i, iCell, k)
- do iVertex = 1, nVertices
- do i=1,grid % vertexDegree
- iCell = cellsOnVertex(i,iVertex)
- if (iCell <= nCells) then
- do k = 1,nVertLevels
- pv_cell(k,iCell) = pv_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
- vorticity_cell(k,iCell) = vorticity_cell(k,iCell) + kiteAreasOnVertex(i, iVertex) * vorticity(k, iVertex) / areaCell(iCell)
- enddo
- endif
- enddo
- enddo
+ !$omp do private(iCell, i, iVertex, j, k)
+ do iCell = 1, nCells
+ do i = 1, grid % nEdgesOnCell % array(iCell)
+ iVertex = grid % verticesOnCell % array(i, iCell)
+ do j = 1, grid % vertexDegree
+ if(cellsOnVertex(j, iVertex) == iCell) then
+ do k = 1, nVertLevels
+ pv_cell(k, iCell) = pv_cell(k, iCell) + kiteAreasOnVertex(j, iVertex) * pv_vertex(k, iVertex) / areaCell(iCell)
+ vorticity_cell(k, iCell) = vorticity_cell(k, iCell) + kiteAreasOnVertex(j, iVertex) * vorticity(k, iVertex) / areaCell(iCell)
+ end do
+ end if
+ end do
+ end do
+ end do
!$omp end do
!
@@ -1771,7 +1343,7 @@
gradPVn(:,:) = 0.0
!$omp end workshare
- !$omp do private(iEdge, k)
+ !$omp do private(iEdge, k)
do iEdge = 1,nEdges
if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
do k = 1,nVertLevels
@@ -1784,7 +1356,7 @@
! Modify PV edge with upstream bias.
!
- !$omp do private(iEdge, k)
+ !$omp do private(iEdge, k)
do iEdge = 1,nEdges
do k = 1,nVertLevels
pv_edge(k,iEdge) = pv_edge(k,iEdge) - config_apvm_upwinding * u(k,iEdge) * dt * gradPVn(k,iEdge)
@@ -1792,6 +1364,8 @@
enddo
!$omp end do
+ !$omp end parallel
+
!
! set pv_edge = fEdge / h_edge at boundary points
!
@@ -1850,7 +1424,7 @@
if(maxval(boundaryEdge).le.0) return
- !$omp do private(iEdge, k)
+ !$omp parallel do private(iEdge, k)
do iEdge = 1,nEdges
do k = 1,nVertLevels
@@ -1860,7 +1434,7 @@
enddo
enddo
- !$omp end do
+ !$omp end parallel do
end subroutine sw_enforce_boundary_edge
</font>
</pre>