<p><b>dwj07@fsu.edu</b> 2011-08-26 11:47:32 -0600 (Fri, 26 Aug 2011)</p><p><br>
Branch Commit<br>
<br>
Splitting out Vmix of tracers into modules<br>
Applying Array masking<br>
<br>
Fusing some loops in module time integration<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-25 22:18:27 UTC (rev 955)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/Makefile 2011-08-26 17:47:32 UTC (rev 956)
@@ -6,6 +6,10 @@
module_OcnHmixTracer.o \
module_OcnHmixTracerDel2.o \
module_OcnHmixTracerDel4.o \
+ module_OcnVmixTracer.o \
+ module_OcnVmixTracerConst.o \
+ module_OcnVmixTracerTanh.o \
+ module_OcnCoriolis.o \
module_time_integration.o \
module_global_diagnostics.o
@@ -18,6 +22,14 @@
module_advection.o:
+module_OcnCoriolis.o:
+
+module_OcnVmixTracerConst.o:
+
+module_OcnVmixTracerTanh.o:
+
+module_OcnVmixTracer.o: module_OcnVmixTracerConst.o module_OcnVmixTracerTanh.o
+
module_OcnHmixTracerDel2.o:
module_OcnHmixTracerDel4.o:
@@ -28,7 +40,7 @@
module_global_diagnostics.o:
-module_mpas_core.o: module_advection.o module_global_diagnostics.o module_test_cases.o module_OcnHmixTracer.o module_OcnHmixTracerDel2.o module_OcnHmixTracerDel4.o module_time_integration.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_time_integration.o
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnCoriolis.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnCoriolis.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnCoriolis.F 2011-08-26 17:47:32 UTC (rev 956)
@@ -0,0 +1,150 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnCoriolis
+!
+!> \brief MPAS ocean horizontal tracer mixing driver
+!> \author Phil Jones
+!> \date 26 July 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the main driver routine for computing
+!> horizontal mixing tendencies. It primarily decides which mixing
+!> parameterizations are being used and calls them individually.
+!> Detailed mixing parameterizations are contained in their own
+!> modules.
+!
+!-----------------------------------------------------------------------
+
+module OcnCoriolis
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnCoriolisTend
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnCoriolisTend
+!
+!> \brief Computes tendency term for horizontal tracer mixing
+!> \author Phil Jones
+!> \date 26 July 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal mixing tendency for tracers
+!> based on current state and user choices of mixing parameterization.
+!> Multiple parameterizations may be chosen and added together. These
+!> tendencies are generally computed by calling the specific routine
+!> for the chosen parameterization, so this routine is primarily a
+!> driver for managing these choices.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnCoriolisTend(tend, pvEdge, hEdge, u, ke, nEdgesSolve, grid)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ hEdge, u, pvEdge, ke
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ integer, intent(in) :: nEdgesSolve
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(inout) :: &
+ tend !< Input/Output: velocity tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iEdge, j, k, cell1, cell2, eoe
+ real (kind=RKIND) :: q, workpv
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnEdge
+ integer, dimension(:), pointer :: maxLeveLEdgeTop, nEdgesOnEdge
+
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
+
+ real (kind=RKIND), dimension(:), pointer :: dcEdge
+
+ cellsOnEdge => grid % cellsOnEdge % array
+ edgesOnEdge => grid % edgesOnEdge % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
+ nEdgesOnEdge => grid % nEdgesOnEdge % array
+ weightsOnEdge => grid % weightsOnEdge % array
+ dcEdge => grid % dcEdge % array
+
+ do iEdge=1,grid % nEdgesSolve
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+
+ do k=1,maxLevelEdgeTop(iEdge)
+
+ q = 0.0
+ do j = 1,nEdgesOnEdge(iEdge)
+ eoe = edgesOnEdge(j,iEdge)
+ workpv = 0.5 * (pvEdge(k,iEdge) + pvEdge(k,eoe))
+ q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * hEdge(k,eoe)
+ end do
+ tend(k,iEdge) = tend(k,iEdge) &
+ + q &
+ - ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
+
+ end do
+ end do
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnCoriolisTend
+
+!***********************************************************************
+
+end module OcnCoriolis
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixTracerDel2.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixTracerDel2.F 2011-08-25 22:18:27 UTC (rev 955)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixTracerDel2.F 2011-08-26 17:47:32 UTC (rev 956)
@@ -119,7 +119,7 @@
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension(:), allocatable:: boundaryMask
+ real (kind=RKIND) :: boundaryMask
real (kind=RKIND), dimension(:,:,:), allocatable:: del2tracer
@@ -156,8 +156,6 @@
!
!-----------------------------------------------------------------
- allocate(boundaryMask(nVertLevels))
-
do iEdge=1,nEdges
!*** get cell pairs that share edge
@@ -167,41 +165,36 @@
invAreaCell1 = 1.d0/areaCell(cell1)
invAreaCell2 = 1.d0/areaCell(cell2)
+ !*** compute \kappa_2 </font>
<font color="red">abla \phi on edge
+
+ do k=1,maxLevelEdgeTop(iEdge)
!*** compute a boundary mask at each vertical level
!*** on this edge
- do k=1,maxLevelEdgeTop(iEdge)
if (boundaryEdge(k,iEdge) == 1) then
- boundaryMask(k) = 0.d0
+ boundaryMask = 0.d0
else
- boundaryMask(k) = 1.d0
+ boundaryMask = 1.d0
endif
- end do
- !*** compute \kappa_2 </font>
<font color="red">abla \phi on edge
+ do iTracer=1,nTracers
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,nTracers
+ ! div(h \kappa_2 </font>
<font color="red">abla \phi) at cell center
- ! div(h \kappa_2 </font>
<font color="gray">abla \phi) at cell center
+ flux = dvEdge(iEdge) * hEdge(k,iEdge) * eddyDiff2 * &
+ ( tracers(iTracer,k,cell2) &
+ - tracers(iTracer,k,cell1))/dcEdge(iEdge)* &
+ boundaryMask
- flux = dvEdge(iEdge) * hEdge(k,iEdge) * eddyDiff2 * &
- ( tracers(iTracer,k,cell2) &
- - tracers(iTracer,k,cell1))/dcEdge(iEdge)* &
- boundaryMask(k)
+ tend(iTracer,k,cell1) = &
+ tend(iTracer,k,cell1) + flux * invAreaCell1
+ tend(iTracer,k,cell2) = &
+ tend(iTracer,k,cell2) - flux * invAreaCell2
- 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
-
end do
- deallocate(boundaryMask)
-
!--------------------------------------------------------------------
end subroutine OcnHMixTracerDel2Tend
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixTracerDel4.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixTracerDel4.F 2011-08-25 22:18:27 UTC (rev 955)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnHmixTracerDel4.F 2011-08-26 17:47:32 UTC (rev 956)
@@ -121,7 +121,7 @@
integer, dimension(:,:), pointer :: cellsOnEdge
- real (kind=RKIND), dimension(:), allocatable:: boundaryMask
+ real (kind=RKIND) :: boundaryMask
real (kind=RKIND), dimension(:,:,:), allocatable:: del2tracer
@@ -160,7 +160,6 @@
!
!-----------------------------------------------------------------
- allocate(boundaryMask(nVertLevels))
allocate(del2tracer(nTracers,nVertLevels,nCells+1))
del2tracer(:,:,:) = 0.d0
@@ -175,34 +174,32 @@
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
+
+
+ do k=1,maxLevelEdgeTop(iEdge)
!*** compute a boundary mask at each vertical level
!*** on this edge
- do k=1,maxLevelEdgeTop(iEdge)
if (boundaryEdge(k,iEdge) == 1) then
- boundaryMask(k) = 0.d0
+ boundaryMask = 0.d0
else
- boundaryMask(k) = 1.d0
+ boundaryMask = 1.d0
endif
- end do
+ do iTracer=1,nTracers
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,nTracers
+ flux = dvEdge(iEdge)*hEdge(k,iEdge) * &
+ (tracers(iTracer,k,cell2) - &
+ tracers(iTracer,k,cell1))/dcEdge(iEdge)* &
+ boundaryMask
- flux = dvEdge(iEdge)*hEdge(k,iEdge) * &
- (tracers(iTracer,k,cell2) - &
- tracers(iTracer,k,cell1))/dcEdge(iEdge)* &
- boundaryMask(k)
+ del2tracer(iTracer,k,cell1) = &
+ del2tracer(iTracer,k,cell1) + flux
+
+ del2tracer(iTracer,k,cell2) = &
+ del2tracer(iTracer,k,cell2) - flux
- del2tracer(iTracer,k,cell1) = &
- del2tracer(iTracer,k,cell1) + flux
-
- del2tracer(iTracer,k,cell2) = &
- del2tracer(iTracer,k,cell2) - flux
-
+ end do
end do
- end do
-
end do
do iCell = 1,nCells
@@ -227,35 +224,33 @@
invAreaCell1 = 1.0 / areaCell(cell1)
invAreaCell2 = 1.0 / areaCell(cell2)
+
+
+ do k=1,maxLevelEdgeTop(iEdge)
!*** compute a boundary mask at each vertical level
!*** on this edge
- do k=1,maxLevelEdgeTop(iEdge)
if (boundaryEdge(k,iEdge) == 1) then
- boundaryMask(k) = 0.d0
+ boundaryMask = 0.d0
else
- boundaryMask(k) = 1.d0
+ boundaryMask = 1.d0
endif
- end do
+ do iTracer=1,nTracers
- do k=1,maxLevelEdgeTop(iEdge)
- do iTracer=1,nTracers
+ flux = dvEdge(iEdge)* eddyDiff4 * &
+ (del2tracer(iTracer,k,cell2) - &
+ del2tracer(iTracer,k,cell1))/dcEdge(iEdge) * &
+ boundaryMask
- flux = dvEdge(iEdge)* eddyDiff4 * &
- (del2tracer(iTracer,k,cell2) - &
- del2tracer(iTracer,k,cell1))/dcEdge(iEdge) * &
- boundaryMask(k)
-
- tend(iTracer,k,cell1) = &
- tend(iTracer,k,cell1) - flux * invAreaCell1
- tend(iTracer,k,cell2) = &
- tend(iTracer,k,cell2) + flux * invAreaCell2
-
+ tend(iTracer,k,cell1) = &
+ tend(iTracer,k,cell1) - flux * invAreaCell1
+ tend(iTracer,k,cell2) = &
+ tend(iTracer,k,cell2) + flux * invAreaCell2
+ enddo
enddo
- enddo
end do
- deallocate(del2tracer, boundaryMask)
+ deallocate(del2tracer)
!--------------------------------------------------------------------
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracer.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracer.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracer.F 2011-08-26 17:47:32 UTC (rev 956)
@@ -0,0 +1,176 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVmixTracer
+!
+!> \brief MPAS ocean horizontal tracer mixing driver
+!> \author Phil Jones
+!> \date 26 July 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the main driver routine for computing
+!> horizontal mixing tendencies. It primarily decides which mixing
+!> parameterizations are being used and calls them individually.
+!> Detailed mixing parameterizations are contained in their own
+!> modules.
+!
+!-----------------------------------------------------------------------
+
+module OcnVmixTracer
+
+ use grid_types
+ use configure
+ use OcnVmixTracerConst
+ use OcnVmixTracerTanh
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnVmixTracerTend, &
+ OcnVmixTracerInit
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVmixTracerTend
+!
+!> \brief Computes tendency term for horizontal tracer mixing
+!> \author Phil Jones
+!> \date 26 July 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal mixing tendency for tracers
+!> based on current state and user choices of mixing parameterization.
+!> Multiple parameterizations may be chosen and added together. These
+!> tendencies are generally computed by calling the specific routine
+!> for the chosen parameterization, so this routine is primarily a
+!> driver for managing these choices.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVmixTracerTend(tend, tracers, hCell, nCellsSolve, grid, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ tracers !< Input: array of tracers
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ hCell !< Input: layer thickness at cell centers
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ integer, intent(in) :: nCellsSolve
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+ tend !< Input/Output: tracer tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(inout) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! call relevant routines for computing tendencies
+ ! note that the user can choose multiple options and the
+ ! tendencies will be added together
+ !
+ !-----------------------------------------------------------------
+
+ call OcnVmixTracerConstTend(tend, tracers, hCell, nCellsSolve, grid, err)
+ call OcnVmixTracerTanhTend(tend, tracers, hCell, nCellsSolve , grid, err)
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixTracerTend
+
+!***********************************************************************
+!
+! routine OcnVmixTracerInit
+!
+!> \brief Initializes ocean tracer horizontal mixing quantities
+!> \author Phil Jones
+!> \date 26 July 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> horizontal tracer mixing in the ocean. Since a variety of
+!> parameterizations are available, this routine primarily calls the
+!> individual init routines for each parameterization.
+!
+!-----------------------------------------------------------------------
+
+
+ subroutine OcnVmixTracerInit(err)
+
+ !--------------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! Output Variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! call individual init routines for each parameterization
+ !
+ !-----------------------------------------------------------------
+
+ call OcnVmixTracerConstInit(err)
+ call OcnVmixTracerTanhInit(err)
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixTracerInit
+
+!***********************************************************************
+
+end module OcnVmixTracer
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracerConst.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracerConst.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracerConst.F 2011-08-26 17:47:32 UTC (rev 956)
@@ -0,0 +1,230 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVmixTracerConst
+!
+!> \brief Ocean horizontal mixing - Laplacian parameterization
+!> \author Phil Jones
+!> \date 1 August 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines for computing horizontal mixing
+!> tendencies using a Laplacian formulation.
+!
+!-----------------------------------------------------------------------
+
+module OcnVmixTracerConst
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnVmixTracerConstTend, &
+ OcnVmixTracerConstInit
+
+ !-------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: &
+ VmixConstOn !< local flag to determine whether Const chosen
+
+ real (kind=RKIND) :: &
+ diff !< base diffusion for vertical mixing
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVmixTracerConstTend
+!
+!> \brief Computes tendency term for Laplacian horizontal tracer mixing
+!> \author Phil Jones
+!> \date 1 August 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal mixing tendency for tracers
+!> based on a Laplacian form for the mixing. This tendency takes the
+!> form </font>
<font color="black">abla \cdot (h \kappa_2 </font>
<font color="gray">abla \phi), where \phi is the
+!> tracer, h is the thickness and \kappa_2 is a base horizontal
+!> diffusivity.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVmixTracerConstTend(tend, tracers, h, nCellsSolve, grid, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ tracers !< Input: array of tracers
+
+ integer, intent(in) :: &
+ nCellsSolve !< Input: number of cells to solve
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h !< Input: thickness at cell centers
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+ tend !< Input/Output: tracer tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! Output Variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: &
+ k, iCell, iTracer, &! loop counters
+ nVertLevels, nTracers ! array sizes
+
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND), dimension(:), allocatable :: vertDiffTop
+
+ real (kind=RKIND), dimension(:,:), allocatable :: fluxVertTop
+
+ !-----------------------------------------------------------------
+ !
+ ! exit if this mixing is not selected
+ !
+ !-----------------------------------------------------------------
+
+ if (.not. VmixConstOn) return
+ err = 0
+
+ !-----------------------------------------------------------------
+ !
+ ! initialize some arrays and sizes
+ !
+ !-----------------------------------------------------------------
+
+ maxLevelCell => grid % maxLevelCell % array
+
+ nTracers = size(tracers,dim=1)
+ nVertLevels = size(tracers, dim=2)
+
+
+ allocate(vertDiffTop(nVertLevels+1))
+ allocate(fluxVertTop(nTracers,nVertLevels+1))
+
+ vertDiffTop = diff
+
+ fluxVertTop(:,1) = 0.0
+
+ do iCell=1,nCellsSolve
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,nTracers
+ ! compute \kappa_v d\phi/dz
+ fluxVertTop(iTracer,k) = vertDiffTop(k) &
+ *(tracers(iTracer,k-1,iCell) - tracers(iTracer,k,iCell))&
+ *2/(h(k-1,iCell) + h(k,iCell))
+ enddo
+ enddo
+
+ fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
+
+ do k=1,maxLevelCell(iCell)
+ do iTracer=1,nTracers
+ ! This is h d/dz(fluxVertTop) but h and
+ ! dz cancel, so reduces to delta(fluxVertTop)
+ tend(iTracer,k,iCell) = tend(iTracer,k,iCell) &
+ +fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
+ enddo
+ enddo
+ enddo
+ ! iCell
+ ! loop
+ deallocate(fluxVertTop, vertDiffTop)
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixTracerConstTend
+
+!***********************************************************************
+!
+! routine OcnVmixTracerConstInit
+!
+!> \brief Initializes ocean tracer Laplacian horizontal mixing
+!> \author Phil Jones
+!> \date 26 July 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> Laplacian horizontal tracer mixing in the ocean.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVmixTracerConstInit(err)
+
+
+ !--------------------------------------------------------------------
+ !
+ ! Output variables
+ !
+ !--------------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !--------------------------------------------------------------------
+ !
+ ! set some local module variables based on input config choices
+ !
+ !--------------------------------------------------------------------
+
+ VmixConstOn = .false.
+ err = 0
+
+ if ( config_vert_diff_type .EQ. 'const' ) then
+ VmixConstOn = .true.
+
+ diff = config_vert_diffusion
+ endif
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixTracerConstInit
+
+!***********************************************************************
+
+end module OcnVmixTracerConst
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Added: branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracerTanh.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracerTanh.F (rev 0)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_OcnVmixTracerTanh.F 2011-08-26 17:47:32 UTC (rev 956)
@@ -0,0 +1,241 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! OcnVmixTracerTanh
+!
+!> \brief Ocean horizontal mixing - Laplacian parameterization
+!> \author Phil Jones
+!> \date 1 August 2011
+!> \version SVN:$Id:$
+!> \details
+!> This module contains routines for computing horizontal mixing
+!> tendencies using a Laplacian formulation.
+!
+!-----------------------------------------------------------------------
+
+module OcnVmixTracerTanh
+
+ use grid_types
+ use configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: OcnVmixTracerTanhTend, &
+ OcnVmixTracerTanhInit
+
+ !-------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: &
+ VmixTanhOn !< local flag to determine whether Tanh chosen
+
+ real (kind=RKIND) :: &
+ diff !< base diffusion for vertical mixing
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine OcnVmixTracerTanhTend
+!
+!> \brief Computes tendency term for Laplacian horizontal tracer mixing
+!> \author Phil Jones
+!> \date 1 August 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the horizontal mixing tendency for tracers
+!> based on a Laplacian form for the mixing. This tendency takes the
+!> form </font>
<font color="black">abla \cdot (h \kappa_2 </font>
<font color="gray">abla \phi), where \phi is the
+!> tracer, h is the thickness and \kappa_2 is a base horizontal
+!> diffusivity.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVmixTracerTanhTend(tend, tracers, hCell, nCellsSolve, grid, err)
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: &
+ tracers !< Input: array of tracers
+
+ integer , intent(in) :: &
+ nCellsSolve !< Input: number of cells
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ hCell !< Input: thickness at cell centers
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+ tend !< Input/Output: tracer tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: &
+ k, iCell, iTracer, &! loop counters
+ nVertLevels, nTracers ! array sizes
+
+ integer, dimension(:), pointer :: maxLevelCell
+
+ real (kind=RKIND), dimension(:), pointer :: zTopZLevel
+
+ real (kind=RKIND), dimension(:), allocatable :: vertDiffTop
+
+ real (kind=RKIND), dimension(:,:), allocatable :: fluxVertTop
+
+
+ !-----------------------------------------------------------------
+ !
+ ! exit if this mixing is not selected
+ !
+ !-----------------------------------------------------------------
+
+ if (.not. VmixTanhOn) return
+ err = 0
+
+ !-----------------------------------------------------------------
+ !
+ ! initialize some arrays and sizes
+ !
+ !-----------------------------------------------------------------
+
+ maxLevelCell => grid % maxLevelCell % array
+ zTopZLevel => grid % zTopZLevel % array
+
+ nTracers = size(tracers,dim=1)
+ nVertLevels = size(tracers, dim=2)
+
+ allocate(vertDiffTop(nVertLevels+1))
+ allocate(fluxVertTop(nTracers,nVertLevels+1))
+
+ do k=1,nVertLevels+1
+ vertDiffTop(k) = -(config_vmixTanhDiffMax-config_vmixTanhDiffMin)/2.0 &
+ *tanh(-(zTopZLevel(k)-config_vmixTanhZMid) &
+ /config_vmixTanhZWidth) &
+ +(config_vmixTanhDiffMax+config_vmixTanhDiffMin)/2
+ enddo
+
+ fluxVertTop(:,1) = 0.0
+ do iCell=1,nCellsSolve
+ do k=2,maxLevelCell(iCell)
+ do iTracer=1,nTracers
+ ! compute \kappa_v
+ ! d\phi/dz
+ fluxVertTop(iTracer,k) = vertDiffTop(k) &
+ *(tracers(iTracer,k-1,iCell) - tracers(iTracer,k,iCell))&
+ *2/(hCell(k-1,iCell) + hCell(k,iCell))
+ enddo
+ enddo
+ fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
+ do k=1,maxLevelCell(iCell)
+ do iTracer=1,nTracers
+ ! This is h d/dz(fluxVertTop) but h and
+ ! dz cancel, so reduces to delta(fluxVertTop)
+ tend(iTracer,k,iCell) = tend(iTracer,k,iCell) &
+ +fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
+ enddo
+ enddo
+ enddo
+ ! iCell
+ ! loop
+ deallocate(fluxVertTop, vertDiffTop)
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixTracerTanhTend
+
+!***********************************************************************
+!
+! routine OcnVmixTracerTanhInit
+!
+!> \brief Initializes ocean tracer Laplacian horizontal mixing
+!> \author Phil Jones
+!> \date 26 July 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine initializes a variety of quantities related to
+!> Laplacian horizontal tracer mixing in the ocean.
+!
+!-----------------------------------------------------------------------
+
+ subroutine OcnVmixTracerTanhInit(err)
+
+ !--------------------------------------------------------------------
+ !
+ ! Output Variables
+ !
+ !--------------------------------------------------------------------
+
+ integer, intent(out) :: err
+
+ !--------------------------------------------------------------------
+ !
+ ! set some local module variables based on input config choices
+ !
+ !--------------------------------------------------------------------
+
+ VmixTanhOn = .false.
+ err = 0
+
+ if ( config_vert_diff_type .EQ. 'tanh' ) then
+ if (config_vert_grid_type.ne.'zlevel') then
+ write(0,*) 'Abort: config_vert_diff_type.eq.tanh may only', &
+ ' use config_vert_grid_type of zlevel at this time'
+ err = 1
+ endif
+
+ VmixTanhOn = .true.
+
+ diff = config_vert_diffusion
+ endif
+
+ !--------------------------------------------------------------------
+
+ end subroutine OcnVmixTracerTanhInit
+
+!***********************************************************************
+
+end module OcnVmixTracerTanh
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Modified: branches/ocean_projects/performance/mpas/src/core_ocean/module_mpas_core.F
===================================================================
--- branches/ocean_projects/performance/mpas/src/core_ocean/module_mpas_core.F 2011-08-25 22:18:27 UTC (rev 955)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_mpas_core.F 2011-08-26 17:47:32 UTC (rev 956)
@@ -2,8 +2,9 @@
use mpas_framework
use dmpar
- use test_cases
+ use test_cases
use OcnHmixTracer
+ use OcnVmixTracer
type (io_output_object) :: restart_obj
integer :: restart_frame
@@ -24,6 +25,7 @@
real (kind=RKIND) :: dt
type (block_type), pointer :: block
type (dm_info) :: dminfo
+ integer :: err
if (.not. config_do_restart) call setup_sw_test_case(domain)
@@ -55,7 +57,12 @@
!***
call OcnHmixTracerInit
+ call OcnVmixTracerInit(err)
+ if(err .ne. 0) then
+ call dmpar_abort(dminfo)
+ endif
+
! mrp 100316 In order for this to work, we need to pass domain % dminfo as an
! input arguement into mpas_init. Ask about that later. For now, there will be
! no initial statistics write.
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-25 22:18:27 UTC (rev 955)
+++ branches/ocean_projects/performance/mpas/src/core_ocean/module_time_integration.F 2011-08-26 17:47:32 UTC (rev 956)
@@ -10,6 +10,8 @@
use spline_interpolation
use timer
use OcnHmixTracer
+ use OcnVmixTracer
+ use OcnCoriolis
contains
@@ -369,15 +371,10 @@
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
- tend_h(k,cell2) = tend_h(k,cell2) + flux
+ 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
- do iCell=1,nCells
- do k=1,nVertLevels
- tend_h(k,iCell) = tend_h(k,iCell) / areaCell(iCell)
- end do
- end do
elseif (config_vert_grid_type.eq.'zlevel') then
@@ -386,24 +383,16 @@
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
- tend_h(k,cell2) = tend_h(k,cell2) + flux
+ 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) / areaCell(iCell)
+ tend_h(1,iCell) = tend_h(1,iCell) + wTop(2,iCell)
end do
-
- endif ! config_vert_grid_type
-
- !
- ! height tendency: vertical advection term -d/dz(hw)
- !
- ! Vertical advection computed for top layer of a z grid only.
- if (config_vert_grid_type.eq.'zlevel') then
- do iCell=1,nCells
- tend_h(1,iCell) = tend_h(1,iCell) + wTop(2,iCell)
- end do
endif ! coordinate type
call timer_stop("height")
@@ -513,6 +502,7 @@
allocate(delsq_vorticity(nVertLevels, nVertices+1))
delsq_u(:,:) = 0.0
+ delsq_circulation(:,:) = 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
@@ -527,27 +517,21 @@
( divergence(k,cell2) - divergence(k,cell1) ) / dcEdge(iEdge) &
-( vorticity(k,vertex2) - vorticity(k,vertex1)) / dvEdge(iEdge)
+ !---- DWJ Testing Comments for loop merging
end do
end do
! vorticity using </font>
<font color="red">abla^2 u
- delsq_circulation(:,:) = 0.0
do iEdge=1,nEdges
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)
+ - dcEdge(iEdge) * delsq_u(k,iEdge) / areaTriangle(vertex1)
delsq_circulation(k,vertex2) = delsq_circulation(k,vertex2) &
- + dcEdge(iEdge) * delsq_u(k,iEdge)
+ + dcEdge(iEdge) * delsq_u(k,iEdge) / areaTriangle(vertex2)
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="gray">abla^2 u
delsq_divergence(:,:) = 0.0
@@ -556,17 +540,11 @@
cell2 = cellsOnEdge(2,iEdge)
do k=1,maxLevelEdgeTop(iEdge)
delsq_divergence(k,cell1) = delsq_divergence(k,cell1) &
- + delsq_u(k,iEdge)*dvEdge(iEdge)
+ + delsq_u(k,iEdge)*dvEdge(iEdge) / areaCell(cell1)
delsq_divergence(k,cell2) = delsq_divergence(k,cell2) &
- - delsq_u(k,iEdge)*dvEdge(iEdge)
+ - delsq_u(k,iEdge)*dvEdge(iEdge) / areaCell(cell2)
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) )
@@ -595,29 +573,13 @@
end if
call timer_stop("vel dissipation")
- call timer_start("coriolis")
!
! velocity tendency: nonlinear Coriolis term and grad of kinetic energy
!
- do iEdge=1,grid % nEdgesSolve
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
+ call timer_start("coriolis")
- do k=1,maxLevelEdgeTop(iEdge)
+ call OcnCoriolisTend(tend_u, pv_edge, h_edge, u, ke, nEdgesSolve, grid)
- q = 0.0
- do j = 1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- workpv = 0.5 * (pv_edge(k,iEdge) + pv_edge(k,eoe))
- q = q + weightsOnEdge(j,iEdge) * u(k,eoe) * workpv * h_edge(k,eoe)
- end do
- tend_u(k,iEdge) = tend_u(k,iEdge) &
- + q &
- - ( ke(k,cell2) - ke(k,cell1) ) / dcEdge(iEdge)
-
- end do
- end do
-
call timer_stop("coriolis")
call timer_start("vel forcing")
!
@@ -757,6 +719,8 @@
real (kind=RKIND), dimension(:,:,:),allocatable :: &
tracersA, tend_trA
+ integer :: err
+
u => s % u % array
h => s % h % array
boundaryCell=> grid % boundaryCell % array
@@ -1119,54 +1083,13 @@
! tracer tendency: vertical diffusion h d/dz( \kappa_v d\phi/dz)
!
call timer_start("scalar vdiff")
- allocate(vertDiffTop(nVertLevels+1))
- if (config_vert_diff_type.eq.'const') then
- vertDiffTop = config_vert_diffusion
- elseif (config_vert_diff_type.eq.'tanh') then
- if (config_vert_grid_type.ne.'zlevel') then
- write(0,*) 'Abort: config_vert_diff_type.eq.tanh may only', &
- ' use config_vert_grid_type of zlevel at this time'
- call dmpar_abort(dminfo)
- endif
-
- do k=1,nVertLevels+1
- vertDiffTop(k) = -(config_vmixTanhDiffMax-config_vmixTanhDiffMin)/2.0 &
- *tanh(-(zTopZLevel(k)-config_vmixTanhZMid) &
- /config_vmixTanhZWidth) &
- + (config_vmixTanhDiffMax+config_vmixTanhDiffMin)/2
- enddo
- else
- write(0,*) 'Abort: unrecognized config_vert_diff_type'
+ call OcnVmixTracerTend(tend_tr, tracers, h, nCellsSolve, grid, err)
+
+ if(err .ne. 0) then
call dmpar_abort(dminfo)
endif
-
- allocate(fluxVertTop(num_tracers,nVertLevels+1))
- fluxVertTop(:,1) = 0.0
- do iCell=1,nCellsSolve
-
- do k=2,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- ! compute \kappa_v d\phi/dz
- fluxVertTop(iTracer,k) = vertDiffTop(k) &
- * (tracers(iTracer,k-1,iCell) - tracers(iTracer,k,iCell) )&
- * 2 / (h(k-1,iCell) + h(k,iCell))
- enddo
- enddo
- fluxVertTop(:,maxLevelCell(iCell)+1) = 0.0
-
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- ! This is h d/dz( fluxVertTop) but h and dz cancel, so
- ! reduces to delta( fluxVertTop)
- tend_tr(iTracer,k,iCell) = tend_tr(iTracer,k,iCell) &
- + fluxVertTop(iTracer,k) - fluxVertTop(iTracer,k+1)
- enddo
- enddo
-
- enddo ! iCell loop
- deallocate(fluxVertTop, vertDiffTop)
-
call timer_stop("scalar vdiff")
+
call timer_stop("scalar")
! print some diagnostics - for debugging
! print *, 'after vertical mixing',&
@@ -1700,18 +1623,11 @@
boundaryEdge => grid % boundaryEdge % array
tend_u => tend % u % array
- if(maxval(boundaryEdge).le.0) return
+! if(maxval(boundaryEdge).le.0) return
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
-
- if(boundaryEdge(k,iEdge).eq.1) then
- tend_u(k,iEdge) = 0.0
- endif
-
- enddo
- enddo
-
+ where(boundaryEdge .eq. 1)
+ tend_u = 0.0
+ end where
end subroutine enforce_boundaryEdge
</font>
</pre>