<p><b>dwj07@fsu.edu</b> 2011-11-04 15:01:33 -0600 (Fri, 04 Nov 2011)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Some performance optimizations.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_mpas_core.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_mpas_core.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -715,8 +715,8 @@
maxLevelVertexTop(nVertices+1) = 0
! set boundary edge
- boundaryEdge=1
- edgeMask=0
+ boundaryEdge(:,1:nEdges+1)=1
+ edgeMask(:,1:nEdges+1)=0
do iEdge=1,nEdges
boundaryEdge(1:maxLevelEdgeTop(iEdge),iEdge)=0
edgeMask(1:maxLevelEdgeTop(iEdge),iEdge)=1
@@ -725,10 +725,10 @@
!
! Find cells and vertices that have an edge on the boundary
!
- boundaryCell = 0
- boundaryVertex = 0
- cellMask = 1
- vertexMask = 1
+ boundaryCell(:,1:nCells+1) = 0
+ cellMask(:,1:nCells+1) = 1
+ boundaryVertex(:,1:nVertices+1) = 0
+ vertexMask(:,1:nVertices+1) = 1
do iEdge=1,nEdges
do k=1,nVertLevels
if (boundaryEdge(k,iEdge).eq.1) then
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hadv2.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hadv2.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hadv2.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -42,7 +42,7 @@
!
!--------------------------------------------------------------------
- logical :: hadv2On !< Flag to turn on/off 2nd order hadv
+ integer :: hadv2On !< Flag to turn on/off 2nd order hadv
!***********************************************************************
@@ -110,7 +110,7 @@
integer, dimension(:), pointer :: maxLevelEdgeTop
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND) :: flux, tracer_edge
+ real (kind=RKIND) :: flux, tracer_edge, invAreaCell1, invAreaCell2, r_tmp
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
@@ -136,12 +136,16 @@
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
do k=1,maxLevelEdgeTop(iEdge)
+ r_tmp = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
do iTracer=1,num_tracers
tracer_edge = 0.5 * (tracers(iTracer,k,cell1) + tracers(iTracer,k,cell2))
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge) * tracer_edge
- tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - flux/areaCell(cell1)
- tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + flux/areaCell(cell2)
+ flux = r_tmp * tracer_edge
+ tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - flux * invAreaCell1
+ tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + flux * invAreaCell2
end do
end do
end do
@@ -176,10 +180,10 @@
integer, intent(out) :: err !< Output: Error flag
err = 0
- hadv2On = .false.
+ hadv2On = 0
if (config_tracer_adv_order == 2) then
- hadv2On = .true.
+ hadv2On = 1
end if
!--------------------------------------------------------------------
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hadv3.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hadv3.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hadv3.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -42,7 +42,7 @@
!
!--------------------------------------------------------------------
- logical :: hadv3On !< Flag to turn on/off 3rd order hadv
+ integer :: hadv3On !< Flag to turn on/off 3rd order hadv
real (kind=RKIND) :: coef_3rd_order !< Coefficient for 3rd order hadv
!***********************************************************************
@@ -223,10 +223,10 @@
integer, intent(out) :: err !< Output: error flag
err = 0
- hadv3On = .false.
+ hadv3On = 0
if (config_tracer_adv_order == 3) then
- hadv3On = .true.
+ hadv3On = 1
coef_3rd_order = 1.0
if (config_monotonic) coef_3rd_order = 0.25
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hadv4.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hadv4.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hadv4.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -42,7 +42,7 @@
!
!--------------------------------------------------------------------
- logical :: hadv4On !< Flag to turning on/off 4th order hadv
+ integer :: hadv4On !< Flag to turning on/off 4th order hadv
!***********************************************************************
@@ -210,10 +210,10 @@
integer, intent(out) :: err !< Output: Error flag
err = 0
- hadv4On = .false.
+ hadv4On = 0
if (config_tracer_adv_order == 4) then
- hadv4On = .true.
+ hadv4On = 1
end if
!--------------------------------------------------------------------
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hmix_del2.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hmix_del2.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -44,7 +44,7 @@
!
!--------------------------------------------------------------------
- logical :: del2On
+ integer :: del2On
real (kind=RKIND) :: eddyDiff2
@@ -113,10 +113,10 @@
integer, dimension(:,:), allocatable :: boundaryMask
integer, dimension(:), pointer :: maxLevelEdgeTop
- integer, dimension(:,:), pointer :: cellsOnEdge, boundaryEdge
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgeMask
real (kind=RKIND) :: invAreaCell1, invAreaCell2
- real (kind=RKIND) :: tracer_turb_flux, flux
+ real (kind=RKIND) :: tracer_turb_flux, flux, r_tmp
real (kind=RKIND), dimension(:), pointer :: areaCell, dvEdge, dcEdge
real (kind=RKIND), dimension(:), pointer :: meshScalingDel2
@@ -139,7 +139,7 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
cellsOnEdge => grid % cellsOnEdge % array
- boundaryEdge => grid % boundaryEdge % array
+ edgeMask => grid % edgeMask % array
areaCell => grid % areaCell % array
dvEdge => grid % dvEdge % array
dcEdge => grid % dcEdge % array
@@ -148,34 +148,28 @@
!
! compute a boundary mask to enforce insulating boundary conditions in the horizontal
!
- allocate(boundaryMask(nVertLevels, nEdges+1))
- boundaryMask = 1.0
- where(boundaryEdge.eq.1) boundaryMask=0.0
-
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
invAreaCell1 = 1.0/areaCell(cell1)
invAreaCell2 = 1.0/areaCell(cell2)
+ r_tmp = meshScalingDel2(iEdge) * eddyDiff2 * dvEdge(iEdge) / dcEdge(iEdge)
+
do k=1,maxLevelEdgeTop(iEdge)
do iTracer=1,num_tracers
! \kappa_2 </font>
<font color="red">abla \phi on edge
- tracer_turb_flux = meshScalingDel2(iEdge) * eddyDiff2 &
- *( tracers(iTracer,k,cell2) &
- - tracers(iTracer,k,cell1))/dcEdge(iEdge)
+ tracer_turb_flux = tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)
! 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)
+ flux = h_edge(k,iEdge) * tracer_turb_flux * edgeMask(k, iEdge) * r_tmp
+
tend(iTracer,k,cell1) = tend(iTracer,k,cell1) + flux * invAreaCell1
tend(iTracer,k,cell2) = tend(iTracer,k,cell2) - flux * invAreaCell2
end do
end do
end do
-
- deallocate(boundaryMask)
!--------------------------------------------------------------------
end subroutine ocn_tracer_hmix_del2_tend!}}}
@@ -208,10 +202,10 @@
err = 0
- del2on = .false.
+ del2on = 0
if ( config_h_tracer_eddy_diff2 > 0.0 ) then
- del2On = .true.
+ del2On = 1
eddyDiff2 = config_h_tracer_eddy_diff2
endif
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hmix_del4.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_hmix_del4.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -44,7 +44,7 @@
!
!--------------------------------------------------------------------
- logical :: Del4On
+ integer :: Del4On
real (kind=RKIND) :: eddyDiff4
@@ -110,12 +110,10 @@
integer :: iEdge, nEdges, num_tracers, nVertLevels, nCells
integer :: iTracer, k, iCell, cell1, cell2
- integer, dimension(:,:), allocatable :: boundaryMask
-
integer, dimension(:), pointer :: maxLevelEdgeTop, maxLevelCell
- integer, dimension(:,:), pointer :: boundaryEdge, cellsOnEdge
+ integer, dimension(:,:), pointer :: edgeMask, cellsOnEdge
- real (kind=RKIND) :: invAreaCell1, invAreaCell2, r, tracer_turb_flux, flux
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2, tracer_turb_flux, flux, r_tmp
real (kind=RKIND), dimension(:,:,:), allocatable :: delsq_tracer
@@ -141,7 +139,7 @@
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
maxLevelCell => grid % maxLevelCell % array
- boundaryEdge => grid % boundaryEdge % array
+ edgeMask => grid % edgeMask % array
cellsOnEdge => grid % cellsOnEdge % array
dcEdge => grid % dcEdge % array
@@ -149,10 +147,6 @@
areaCell => grid % areaCell % array
meshScalingDel4 => grid % meshScalingDel4 % array
- allocate(boundaryMask(nVertLevels, nEdges+1))
- boundaryMask = 1.0
- where(boundaryEdge.eq.1) boundaryMask=0.0
-
allocate(delsq_tracer(num_tracers,nVertLevels, nCells+1))
delsq_tracer(:,:,:) = 0.
@@ -162,29 +156,24 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+ r_tmp = dvEdge(iEdge) / dcEdge(iEdge)
+
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
+
do k=1,maxLevelEdgeTop(iEdge)
do iTracer=1,num_tracers
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)
+ + edgeMask(k,iEdge) * r_tmp * h_edge(k,iEdge) &
+ *(tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) * invAreaCell1
+
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)
+ - edgeMask(k,iEdge) * r_tmp * h_edge(k,iEdge) &
+ *(tracers(iTracer,k,cell2) - tracers(iTracer,k,cell1)) * invAreaCell2
end do
end do
end do
- do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- delsq_tracer(iTracer,k,iCell) = delsq_tracer(iTracer,k,iCell) * r
- end do
- end do
- end do
-
! second del2: div(h </font>
<font color="gray">abla [delsq_tracer]) at cell center
do iEdge=1,grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
@@ -192,17 +181,15 @@
invAreaCell1 = 1.0 / areaCell(cell1)
invAreaCell2 = 1.0 / areaCell(cell2)
+ r_tmp = meshScalingDel4(iEdge) * eddyDiff4 * dvEdge(iEdge) / dcEdge(iEdge)
+
do k=1,maxLevelEdgeTop(iEdge)
do iTracer=1,num_tracers
- tracer_turb_flux = meshScalingDel4(iEdge) * eddyDiff4 &
- *( delsq_tracer(iTracer,k,cell2) &
- - delsq_tracer(iTracer,k,cell1))/dcEdge(iEdge)
- flux = dvEdge (iEdge) * tracer_turb_flux
+ tracer_turb_flux = delsq_tracer(iTracer,k,cell2) - delsq_tracer(iTracer,k,cell1)
+ flux = tracer_turb_flux * edgeMask(k,iEdge) * r_tmp
- tend(iTracer,k,cell1) = tend(iTracer,k,cell1) &
- - flux * invAreaCell1 * boundaryMask(k,iEdge)
- tend(iTracer,k,cell2) = tend(iTracer,k,cell2) &
- + flux * invAreaCell2 * boundaryMask(k,iEdge)
+ tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - flux * invAreaCell1
+ tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + flux * invAreaCell2
enddo
enddo
@@ -240,10 +227,10 @@
integer, intent(out) :: err !< Output: error flag
err = 0
- Del4on = .false.
+ Del4on = 0
if ( config_h_tracer_eddy_diff4 > 0.0 ) then
- Del4On = .true.
+ Del4On = 1
eddyDiff4 = config_h_tracer_eddy_diff4
endif
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_spline2.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_spline2.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_spline2.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -16,7 +16,6 @@
use mpas_grid_types
use mpas_configure
- use mpas_timer
implicit none
private
@@ -43,7 +42,7 @@
!
!--------------------------------------------------------------------
- logical :: spline2On
+ integer :: spline2On
!***********************************************************************
@@ -125,8 +124,6 @@
if(.not.spline2On) return
! Compute tracerTop using linear interpolation.
- call mpas_timer_start("compute_scalar_tend-vert adv spline 2")
-
nCells = grid % nCells
nCellsSolve = grid % nCellsSolve
nVertLevels = grid % nVertLevels
@@ -162,7 +159,6 @@
deallocate(tracerTop)
- call mpas_timer_stop("compute_scalar_tend-vert adv spline 2")
!--------------------------------------------------------------------
end subroutine ocn_tracer_vadv_spline2_tend!}}}
@@ -195,10 +191,10 @@
err = 0
- spline2On = .false.
+ spline2On = 0
if(config_vert_tracer_adv_order.eq.2) then
- spline2On = .true.
+ spline2On = 1
endif
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_spline3.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_spline3.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_spline3.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -16,7 +16,6 @@
use mpas_grid_types
use mpas_configure
- use mpas_timer
use mpas_spline_interpolation
implicit none
@@ -44,7 +43,7 @@
!
!--------------------------------------------------------------------
- logical :: spline3On
+ integer :: spline3On
!***********************************************************************
@@ -129,8 +128,6 @@
if(.not.spline3On) return
! Compute tracerTop using linear interpolation.
- call mpas_timer_start("compute_scalar_tend-vert adv spline 3")
-
nCells = grid % nCells
nCellsSolve = grid % nCellsSolve
nVertLevels = grid % nVertLevels
@@ -190,8 +187,6 @@
deallocate(tracer2ndDer)
deallocate(tracersIn,tracersOut, posZMidZLevel, posZTopZLevel)
deallocate(tracerTop)
-
- call mpas_timer_stop("compute_scalar_tend-vert adv spline 3")
!--------------------------------------------------------------------
end subroutine ocn_tracer_vadv_spline3_tend!}}}
@@ -224,10 +219,10 @@
err = 0
- spline3On = .false.
+ spline3On = 0
if(config_vert_tracer_adv_order.eq.3) then
- spline3On = .true.
+ spline3On = 1
endif
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_stencil2.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_stencil2.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_stencil2.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -16,7 +16,6 @@
use mpas_grid_types
use mpas_configure
- use mpas_timer
implicit none
private
@@ -43,7 +42,7 @@
!
!--------------------------------------------------------------------
- logical :: stencil2On
+ integer :: stencil2On
!***********************************************************************
@@ -124,9 +123,6 @@
if(.not. stencil2On) return
-
- call mpas_timer_start("compute_scalar_tend-vert adv stencil 2")
-
nCells = grid % nCells
nCellsSolve = grid % nCellsSolve
num_tracers = size(tracers, 1)
@@ -142,7 +138,7 @@
do iTracer=1,num_tracers
tracerTop(iTracer,k,iCell) = &
( tracers(iTracer,k-1,iCell) &
- +tracers(iTracer,k ,iCell))/2.0
+ +tracers(iTracer,k ,iCell)) * 0.5
end do
end do
end do
@@ -158,7 +154,6 @@
end do
deallocate(tracerTop)
- call mpas_timer_stop("compute_scalar_tend-vert adv stencil 2")
!--------------------------------------------------------------------
@@ -193,10 +188,10 @@
integer :: err1, err2, err3
err = 0
- stencil2On = .false.
+ stencil2On = 0
if(config_vert_tracer_adv_order.eq.2) then
- stencil2On = .true.
+ stencil2On = 1
endif
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_stencil3.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_stencil3.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_stencil3.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -43,7 +43,7 @@
!
!--------------------------------------------------------------------
- logical :: stencil3On
+ integer :: stencil3On
!***********************************************************************
@@ -134,8 +134,6 @@
hRatioZLevelK => grid % hRatioZLevelK % array
hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
- call mpas_timer_start("compute_scalar_tend-vert adv stencil 3")
-
allocate(tracerTop(num_tracers,nVertLevels+1,nCells))
! Compute tracerTop using 3rd order stencil. This is the same
@@ -181,8 +179,6 @@
end do
deallocate(tracerTop)
- call mpas_timer_stop("compute_scalar_tend-vert adv stencil 3")
-
!--------------------------------------------------------------------
end subroutine ocn_tracer_vadv_stencil3_tend!}}}
@@ -214,10 +210,10 @@
integer, intent(out) :: err !< Output: error flag
err = 0
- stencil3On = .false.
+ stencil3On = 0
if(config_vert_tracer_adv_order.eq.3) then
- stencil3On = .true.
+ stencil3On = 1
endif
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_stencil4.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_stencil4.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_tracer_vadv_stencil4.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -16,7 +16,6 @@
use mpas_grid_types
use mpas_configure
- use mpas_timer
implicit none
private
@@ -43,7 +42,7 @@
!
!--------------------------------------------------------------------
- logical :: stencil4On
+ integer :: stencil4On
!***********************************************************************
@@ -126,8 +125,6 @@
if(.not. Stencil4On) return
- call mpas_timer_start("compute_scalar_tend-vert adv stencil 4")
-
nCells = grid % nCells
nCellsSolve = grid % nCellsSolve
num_tracers = size(tracers, 1)
@@ -176,8 +173,6 @@
end do
deallocate(tracerTop)
- call mpas_timer_stop("compute_scalar_tend-vert adv stencil 4")
-
!--------------------------------------------------------------------
end subroutine ocn_tracer_vadv_stencil4_tend!}}}
@@ -209,10 +204,10 @@
integer, intent(out) :: err !< Output: error flag
err = 0
- stencil4On = .false.
+ stencil4On = 0
if(config_vert_tracer_adv_order.eq.4) then
- stencil4On = .true.
+ stencil4On = 1
endif
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -129,22 +129,14 @@
do iEdge=1,grid % nEdgesSolve
- k = maxLevelEdgeTop(iEdge)
+ k = max(maxLevelEdgeTop(iEdge), 1)
- ! efficiency note: it would be nice to avoid this
- ! if within a do. This could be done with
- ! k = max(maxLevelEdgeTop(iEdge),1)
- ! and then tend_u(1,iEdge) is just not used for land cells.
+ ! bottom drag is the same as POP:
+ ! -c |u| u where c is unitless and 1.0e-3.
+ ! see POP Reference guide, section 3.4.4.
- if (k>0) then
- ! bottom drag is the same as POP:
- ! -c |u| u where c is unitless and 1.0e-3.
- ! see POP Reference guide, section 3.4.4.
+ tend(k,iEdge) = tend(k,iEdge)-edgeMask(k,iEdge)*(bottomDragCoef*u(k,iEdge)*sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge))
- tend(k,iEdge) = tend(k,iEdge)-edgeMask(k,iEdge)*(bottomDragCoef*u(k,iEdge)*sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge))
-
- endif
-
enddo
Modified: branches/ocean_projects/performance/src/core_ocean/mpas_ocn_vel_hmix_del4.F
===================================================================
--- branches/ocean_projects/performance/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2011-11-04 18:53:55 UTC (rev 1171)
+++ branches/ocean_projects/performance/src/core_ocean/mpas_ocn_vel_hmix_del4.F 2011-11-04 21:01:33 UTC (rev 1172)
@@ -122,7 +122,8 @@
integer, dimension(:,:), pointer :: cellsOnEdge, verticesOnEdge, edgeMask
- real (kind=RKIND) :: u_diffusion, r
+ real (kind=RKIND) :: u_diffusion, invAreaCell1, invAreaCell2, invAreaTri1, &
+ invAreaTri2, invDcEdge, invDvEdge, r_tmp
real (kind=RKIND), dimension(:), pointer :: dcEdge, dvEdge, areaTriangle, &
meshScalingDel4, areaCell
@@ -155,6 +156,9 @@
allocate(delsq_vorticity(nVertLevels, nVertices+1))
delsq_u(:,:) = 0.0
+ delsq_circulation(:,:) = 0.0
+ delsq_vorticity(:,:) = 0.0
+ delsq_divergence(:,:) = 0.0
! Compute </font>
<font color="black">abla^2 u = </font>
<font color="black">abla divergence + k \times </font>
<font color="gray">abla vorticity
do iEdge=1,grid % nEdges
@@ -163,53 +167,48 @@
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
+ invDcEdge = 1.0 / dcEdge(iEdge)
+ invDvEdge = 1.0 / dvEdge(iEdge)
+
do k=1,maxLevelEdgeTop(iEdge)
delsq_u(k,iEdge) = &
- ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -viscVortCoef &
- *( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+ ( divergence(k,cell2) - divergence(k,cell1) ) * invDcEdge &
+ -viscVortCoef *( vorticity(k,vertex2) - vorticity(k,vertex1)) * invDvEdge
end do
end do
- ! vorticity using </font>
<font color="red">abla^2 u
- delsq_circulation(:,:) = 0.0
do iEdge=1,nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
- do k=1,maxLevelEdgeTop(iEdge)
- delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) &
- - dcEdge(iEdge) * delsq_u(k,iEdge)
- delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
- + dcEdge(iEdge) * delsq_u(k,iEdge)
- end do
- end do
- do iVertex=1,nVertices
- r = 1.0 / areaTriangle(iVertex)
- do k=1,maxLevelVertexBot(iVertex)
- delsq_vorticity(k,iVertex) = delsq_circulation(k,iVertex) * r
- end do
- end do
- ! Divergence using </font>
<font color="red">abla^2 u
- delsq_divergence(:,:) = 0.0
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ invAreaTri1 = 1.0 / areaTriangle(vertex1)
+ invAreaTri2 = 1.0 / areaTriangle(vertex2)
+
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
+
do k=1,maxLevelEdgeTop(iEdge)
- delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
- + delsq_u(k,iEdge)*dvEdge(iEdge)
- delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
- - delsq_u(k,iEdge)*dvEdge(iEdge)
+ ! vorticity using </font>
<font color="blue">abla^2 u
+ r_tmp = dcEdge(iEdge) * delsq_u(k,iEdge)
+
+ delsq_circulation(k,vertex1) = delsq_circulation(k,vertex1) - r_tmp
+ delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) + r_tmp
+
+ delsq_vorticity(k,vertex1) = delsq_vorticity(k,vertex1) - r_tmp * invAreaTri1
+ delsq_vorticity(k,vertex2) = delsq_vorticity(k,vertex2) + r_tmp * invAreaTri2
+
+ ! Divergence using </font>
<font color="red">abla^2 u
+ r_tmp = dvEdge(iEdge) * delsq_u(k,iEdge)
+ delsq_divergence(k, cell1) = delsq_divergence(k,cell1) + r_tmp * invAreaCell1
+ delsq_divergence(k, cell2) = delsq_divergence(k,cell2) - r_tmp * invAreaCell2
+
end do
end do
- do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
- do k = 1,maxLevelCell(iCell)
- delsq_divergence(k,iCell) = delsq_divergence(k,iCell) * r
- end do
- 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) )
@@ -219,20 +218,19 @@
vertex1 = verticesOnEdge(1,iEdge)
vertex2 = verticesOnEdge(2,iEdge)
+ invDcEdge = 1.0 / dcEdge(iEdge)
+ invDvEdge = 1.0 / dvEdge(iEdge)
+ r_tmp = meshScalingDel4(iEdge) * eddyVisc4
+
do k=1,maxLevelEdgeTop(iEdge)
delsq_u(k,iEdge) = &
- ( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
- -( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+ ( divergence(k,cell2) - divergence(k,cell1) ) * invDcEdge &
+ -( vorticity(k,vertex2) - vorticity(k,vertex1)) * invDvEdge
- u_diffusion = ( delsq_divergence(k,cell2) &
- - delsq_divergence(k,cell1) ) / dcEdge(iEdge) &
- -viscVortCoef &
- *( delsq_vorticity(k,vertex2) &
- - delsq_vorticity(k,vertex1) ) / dvEdge(iEdge)
+ u_diffusion = (delsq_divergence(k,cell2) - delsq_divergence(k,cell1)) * invDcEdge &
+ -viscVortCoef * (delsq_vorticity(k,vertex2) - delsq_vorticity(k,vertex1) ) * invDvEdge
- u_diffusion = meshScalingDel4(iEdge) * eddyVisc4 * u_diffusion
-
- tend(k,iEdge) = tend(k,iEdge) - edgeMask(k, iEdge) * u_diffusion
+ tend(k,iEdge) = tend(k,iEdge) - edgeMask(k, iEdge) * u_diffusion * r_tmp
end do
end do
</font>
</pre>