<p><b>dwj07@fsu.edu</b> 2011-08-29 16:03:13 -0600 (Mon, 29 Aug 2011)</p><p><br>
Separating out thickness tendency and velocity forcing tendencies into modules, and applying array masking.<br>
<br>
changing a lot of divisions from module_time_integration to multiplications, to help performance.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/performance/mpas/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/Makefile 2011-08-29 16:48:20 UTC (rev 958)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/Makefile 2011-08-29 22:03:13 UTC (rev 959)
@@ -13,6 +13,8 @@
module_OcnVmixTracerConst.o \
module_OcnVmixTracerTanh.o \
module_OcnCoriolis.o \
+ module_OcnThickness.o \
+ module_OcnVelocityForcing.o \
module_time_integration.o \
module_global_diagnostics.o
@@ -23,6 +25,10 @@
module_test_cases.o:
+module_OcnVelocityForcing.o:
+
+module_OcnThickness.o:
+
module_advection.o:
module_OcnCoriolis.o:
@@ -49,7 +55,7 @@
module_global_diagnostics.o:
-module_mpas_core.o: module_advection.o module_global_diagnostics.o module_test_cases.o module_OcnCoriolis.o module_OcnVmixTracer.o module_OcnVmixTracerConst.o module_OcnVmixTracerTanh.o module_OcnHmixTracer.o module_OcnHmixTracerDel2.o module_OcnHmixTracerDel4.o module_OcnHmixVel.o module_OcnHmixVelDel2.o module_OcnHmixVelDel4.o module_time_integration.o
+module_mpas_core.o: module_advection.o module_global_diagnostics.o module_test_cases.o module_OcnVelocityForcing.o module_OcnThickness.o module_OcnCoriolis.o module_OcnVmixTracer.o module_OcnVmixTracerConst.o module_OcnVmixTracerTanh.o module_OcnHmixTracer.o module_OcnHmixTracerDel2.o module_OcnHmixTracerDel4.o module_OcnHmixVel.o module_OcnHmixVelDel2.o module_OcnHmixVelDel4.o module_time_integration.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnThickness.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnThickness.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnThickness.F 2011-08-29 22:03:13 UTC (rev 959)
@@ -0,0 +1,168 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnThickness
+!
+!> \brief MPAS ocean thickness tendency
+!> \author Doug Jacobsen
+!> \date 29 July 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the main driver routine for computing
+!> the thickness tendencices
+!
+!-----------------------------------------------------------------------
+
+module OcnThickness
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnThicknessTend
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnThicknessTend
+!
+!> \brief Computes tendency term for horizontal tracer mixing
+!> \author Doug Jacobsen
+!> \date 26 July 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the thickness tendency based on current state.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnThicknessTend(tend, hEdge, u, wTop, grid)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ u, wTop, hEdge
+
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ tend !< Input/Output: velocity tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, nEdges, k, cell1, cell2, iCell, nCells, nVertLevels
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2, flux
+
+ integer, dimension(:,:), pointer :: cellsOnEdge
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
+
+ cellsOnEdge => grid % cellsOnEdge % array
+ areaCell => grid % areaCell % array
+ dvEdge => grid % dvEdge % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+
+ nEdges = grid % nEdges
+ nCells = grid % nCells
+ nVertLevels = grid % nVertLevels
+
+ !
+ ! height tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( hu)
+ !
+ ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
+ ! for explanation of divergence operator.
+ !
+ ! for z-level, only compute height tendency for top layer.
+
+ if (config_vert_grid_type.eq.'isopycnal') then
+
+ 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,nVertLevels
+ flux = u(k,iEdge) * dvEdge(iEdge) * hEdge(k,iEdge)
+ tend(k,cell1) = tend(k,cell1) - flux * invAreaCell1
+ tend(k,cell2) = tend(k,cell2) + flux * invAreaCell2
+ end do
+ end do
+
+ elseif (config_vert_grid_type.eq.'zlevel') then
+
+ 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,min(1,maxLevelEdgeTop(iEdge))
+ flux = u(k,iEdge) * dvEdge(iEdge) * hEdge(k,iEdge)
+ tend(k,cell1) = tend(k,cell1) - flux * invAreaCell1
+ tend(k,cell2) = tend(k,cell2) + flux * invAreaCell2
+ end do
+ end do
+ ! height tendency: vertical advection term -d/dz(hw)
+ !
+ ! Vertical advection computed for top layer of a z grid only.
+ do iCell=1,nCells
+ tend(1,iCell) = tend(1,iCell) + wTop(2,iCell)
+ end do
+ endif ! coordinate type
+
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnThicknessTend
+
+!***********************************************************************
+
+end module OcnThickness
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVelocityForcing.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVelocityForcing.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVelocityForcing.F 2011-08-29 22:03:13 UTC (rev 959)
@@ -0,0 +1,152 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVelocityForcing
+!
+!> \brief MPAS ocean velocity forcing
+!> \author Doug Jacobsen
+!> \date 29 July 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the routine for computing
+!> velocity tendencies from forcings. Including bottom drag and
+!> wind stress.
+!>
+!
+!-----------------------------------------------------------------------
+
+module OcnVelocityForcing
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnVelocityForcingTend
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVelocityForcingTend
+!
+!> \brief Computes tendency term velocity based on forcings
+!> \author Doug Jacobsen
+!> \date 26 July 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the tendency for velocity
+!> based on current state, and forcings that are based on wind stress
+!> (from an initial condition), or bottom drag.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVelocityForcingTend(tend, keEdge, hEdge, u, grid)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ hEdge, u, keEdge
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ tend !< Input/Output: velocity tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, k
+ integer, dimension(:), pointer :: maxLevelEdgeTop
+
+ real (kind=RKIND), dimension(:,:), pointer :: u_src
+
+ real (kind=RKIND), parameter :: rho_ref = 1000.0
+
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ u_src => grid % u_src % array
+
+ !
+ ! velocity tendency: forcing and bottom drag
+ !
+ ! mrp 101115 note: in order to include flux boundary conditions, we will need to
+ ! know the bottom edge with nonzero velocity and place the drag there.
+
+ do iEdge=1,grid % nEdgesSolve
+
+ k = maxLevelEdgeTop(iEdge)
+
+ ! 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.
+
+ if (k>0) then
+
+ ! forcing in top layer only
+ tend(1,iEdge) = tend(1,iEdge) &
+ + u_src(1,iEdge)/rho_ref/hEdge(1,iEdge)
+
+ ! 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) &
+ - 1.0e-3*u(k,iEdge) &
+ *sqrt(2.0*keEdge(k,iEdge))/hEdge(k,iEdge)
+
+ endif
+
+ enddo
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVelocityForcingTend
+
+!***********************************************************************
+
+end module OcnVelocityForcing
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_time_integration.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_time_integration.F 2011-08-29 16:48:20 UTC (rev 958)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_time_integration.F 2011-08-29 22:03:13 UTC (rev 959)
@@ -9,10 +9,12 @@
use vector_reconstruction
use spline_interpolation
use timer
- use OcnHmixVel
use OcnHmixTracer
+ use OcnHmixVel
use OcnVmixTracer
use OcnCoriolis
+ use OcnThickness
+ use OcnVelocityForcing
contains
@@ -299,6 +301,9 @@
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2
+ real (kind=RKIND) :: invAreaTri1, invAreaTri2
+ real (kind=RKIND) :: invDCEdge, invDVEdge
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
@@ -357,45 +362,9 @@
!
tend_h = 0.0
- !
- ! height tendency: horizontal advection term -</font>
<font color="gray">abla\cdot ( hu)
- !
- ! See Ringler et al. (2010) jcp paper, eqn 19, 21, and fig. 3.
- ! for explanation of divergence operator.
- !
- ! for z-level, only compute height tendency for top layer.
+ call OcnThicknessTend(tend_h, h_edge, u, wTop, grid)
- if (config_vert_grid_type.eq.'isopycnal') then
-
- 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)
- tend_h(k,cell1) = tend_h(k,cell1) - flux / areaCell(cell1)
- tend_h(k,cell2) = tend_h(k,cell2) + flux / areaCell(cell2)
- end do
- end do
- elseif (config_vert_grid_type.eq.'zlevel') then
-
- do iEdge=1,nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k=1,min(1,maxLevelEdgeTop(iEdge))
- flux = u(k,iEdge) * dvEdge(iEdge) * h_edge(k,iEdge)
- tend_h(k,cell1) = tend_h(k,cell1) - flux / areaCell(cell1)
- tend_h(k,cell2) = tend_h(k,cell2) + flux / areaCell(cell2)
- end do
- end do
- ! height tendency: vertical advection term -d/dz(hw)
- !
- ! Vertical advection computed for top layer of a z grid only.
- do iCell=1,nCells
- tend_h(1,iCell) = tend_h(1,iCell) + wTop(2,iCell)
- end do
- endif ! coordinate type
-
call timer_stop("height")
call timer_start("velocity")
!
@@ -442,29 +411,35 @@
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+
+ invDCEdge = 1.0 / dcEdge(iEdge)
+
do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- - (MontPot(k,cell2) - MontPot(k,cell1))/dcEdge(iEdge)
+ - (MontPot(k,cell2) - MontPot(k,cell1))*invDCEdge
end do
enddo
elseif (config_vert_grid_type.eq.'zlevel') then
do iEdge=1,grid % nEdgesSolve
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+
+ invDCEdge = 1.0 / dcEdge(iEdge)
+
do k=1,maxLevelEdgeTop(iEdge)
tend_u(k,iEdge) = tend_u(k,iEdge) &
- rho0Inv*( pressure(k,cell2) &
- - pressure(k,cell1) )/dcEdge(iEdge)
+ - pressure(k,cell1) )*invDCEdge
end do
enddo
endif
call timer_stop("press grad")
-
call timer_start("vel dissipation")
+
call OcnHmixVelTend(tend_u, divergence, vorticity, grid)
- call timer_stop("vel dissipation")
+ call timer_stop("vel dissipation")
!
! velocity tendency: nonlinear Coriolis term and grad of kinetic energy
!
@@ -474,39 +449,9 @@
call timer_stop("coriolis")
call timer_start("vel forcing")
- !
- ! velocity tendency: forcing and bottom drag
- !
- ! mrp 101115 note: in order to include flux boundary conditions, we will need to
- ! know the bottom edge with nonzero velocity and place the drag there.
- do iEdge=1,grid % nEdgesSolve
+ call OcnVelocityForcingTend(tend_u, ke_edge, h_edge, u, grid)
- k = maxLevelEdgeTop(iEdge)
-
- ! 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.
-
- if (k>0) then
-
- ! forcing in top layer only
- tend_u(1,iEdge) = tend_u(1,iEdge) &
- + u_src(1,iEdge)/rho_ref/h_edge(1,iEdge)
-
- ! 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_u(k,iEdge) = tend_u(k,iEdge) &
- - 1.0e-3*u(k,iEdge) &
- *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
-
- endif
-
- enddo
-
call timer_stop("vel forcing")
call timer_start("vel vmix")
!
@@ -667,12 +612,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)
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_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux*invAreaCell1
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux*invAreaCell2
end do
end do
end do
@@ -683,6 +632,9 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
+
do k=1,maxLevelEdgeTop(iEdge)
d2fdx2_cell1 = 0.0
@@ -725,8 +677,8 @@
end if
!-- update tendency
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux*invAreaCell1
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux*invAreaCell2
enddo
end do
end do
@@ -737,6 +689,9 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
+ invAreaCell1 = 1.0 / areaCell(cell1)
+ invAreaCell2 = 1.0 / areaCell(cell2)
+
do k=1,maxLevelEdgeTop(iEdge)
d2fdx2_cell1 = 0.0
@@ -769,8 +724,8 @@
-(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. )
!-- update tendency
- tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux/areaCell(cell1)
- tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux/areaCell(cell2)
+ tend_tr(iTracer,k,cell1) = tend_tr(iTracer,k,cell1) - flux*invAreaCell1
+ tend_tr(iTracer,k,cell2) = tend_tr(iTracer,k,cell2) + flux*invAreaCell2
enddo
end do
end do
@@ -1043,6 +998,9 @@
real (kind=RKIND) :: coef_3rd_order
real (kind=RKIND) :: r, h1, h2
+ real (kind=RKIND) :: invAreaCell1, invAreaCell2
+ real (kind=RKIND) :: invAreaTri1, invAreaTri2
+ real (kind=RKIND) :: invDCEdge, invDVEdge
h => s % h % array
u => s % u % array
@@ -1229,8 +1187,11 @@
end do
end do
do iVertex=1,nVertices
+
+ invAreaTri1 = 1.0 / areaTriangle(iVertex)
+
do k=1,maxLevelVertexBot(iVertex)
- vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
+ vorticity(k,iVertex) = circulation(k,iVertex) * invAreaTri1
end do
end do
@@ -1241,17 +1202,15 @@
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,maxLevelEdgeBot(iEdge)
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
- divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+ divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)*invAreaCell1
+ divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)*invAreaCell2
enddo
end do
- do iCell = 1,nCells
- r = 1.0 / areaCell(iCell)
- do k = 1,maxLevelCell(iCell)
- divergence(k,iCell) = divergence(k,iCell) * r
- enddo
- enddo
!
! Compute kinetic energy in each cell
@@ -1260,16 +1219,15 @@
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,maxLevelEdgeBot(iEdge)
- ke(k,cell1) = ke(k,cell1) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
- ke(k,cell2) = ke(k,cell2) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0
+ ke(k,cell1) = ke(k,cell1) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0 * invAreaCell1
+ ke(k,cell2) = ke(k,cell2) + 0.25 * dcEdge(iEdge) * dvEdge(iEdge) * u(k,iEdge)**2.0 * invAreaCell2
enddo
end do
- do iCell = 1,nCells
- do k = 1,maxLevelCell(iCell)
- ke(k,iCell) = ke(k,iCell) / areaCell(iCell)
- enddo
- enddo
!
! Compute v (tangential) velocities
@@ -1305,12 +1263,14 @@
! ( this computes pv_vertex at all vertices bounding real cells and distance-1 ghost cells )
!
do iVertex = 1,nVertices
+
+ invAreaTri1 = 1.0 / areaTriangle(iVertex)
do k=1,maxLevelVertexBot(iVertex)
h_vertex = 0.0
do i=1,vertexDegree
h_vertex = h_vertex + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
end do
- h_vertex = h_vertex / areaTriangle(iVertex)
+ h_vertex = h_vertex * invAreaTri1
pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex
end do
@@ -1324,10 +1284,12 @@
do iVertex = 1,nVertices
do i=1,vertexDegree
iCell = cellsOnVertex(i,iVertex)
+
+ invAreaCell1 = 1.0/areaCell(iCell)
do k = 1,maxLevelCell(iCell)
pv_cell(k,iCell) = pv_cell(k,iCell) &
+ kiteAreasOnVertex(i, iVertex) * pv_vertex(k, iVertex) &
- / areaCell(iCell)
+ * invAreaCell1
enddo
enddo
enddo
@@ -1352,10 +1314,13 @@
!
gradPVn(:,:) = 0.0
do iEdge = 1,nEdges
+
+ invDCEdge = 1.0 / dcEdge(iEdge)
+
do k=1,maxLevelEdgeTop(iEdge)
gradPVn(k,iEdge) = ( pv_cell(k,cellsOnEdge(2,iEdge)) &
- pv_cell(k,cellsOnEdge(1,iEdge))) &
- / dcEdge(iEdge)
+ * invDCEdge
enddo
enddo
@@ -1364,10 +1329,12 @@
! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
!
do iEdge = 1,nEdges
+
+ invDVEdge = 1.0 / dvEdge(iEdge)
do k = 1,maxLevelEdgeBot(iEdge)
gradPVt(k,iEdge) = ( pv_vertex(k,verticesOnEdge(2,iEdge)) &
- pv_vertex(k,verticesOnEdge(1,iEdge))) &
- /dvEdge(iEdge)
+ * invDVEdge
enddo
enddo
@@ -1475,9 +1442,12 @@
! bottom is zero.
wTop(1,iCell) = 0.0
wTop(maxLevelCell(iCell)+1,iCell) = 0.0
+
+ invAreaCell1 = 1.0/areaCell(iCell)
+
do k=maxLevelCell(iCell),2,-1
wTop(k,iCell) = wTop(k+1,iCell) &
- - div_u(k,iCell)/areaCell(iCell)*h(k,iCell)
+ - div_u(k,iCell)*h(k,iCell) * invAreaCell1
end do
end do
deallocate(div_u)
</font>
</pre>