<p><b>ringler@lanl.gov</b> 2011-04-27 11:40:33 -0600 (Wed, 27 Apr 2011)</p><p><br>
removed extra (not relevant) computes in compute_solve_diagnostics<br>
removed unused arrays from Registry<br>
added user-specified nCheck and nQuit parameters in elliptic solves<br>
added profiling to module_time_integration<br>
</p><hr noshade><pre><font color="gray">Modified: branches/pv_based_swm/mpas/src/core_pvsw/Registry
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-04-27 16:10:24 UTC (rev 803)
+++ branches/pv_based_swm/mpas/src/core_pvsw/Registry 2011-04-27 17:40:33 UTC (rev 804)
@@ -145,22 +145,10 @@
# Diagnostic fields: only written to output
var persistent real v ( nVertLevels nEdges Time ) 2 o v state - -
var persistent real divergence ( nVertLevels nCells Time ) 2 o divergence state - -
-var persistent real vorticity ( nVertLevels nVertices Time ) 2 o vorticity state - -
-var persistent real vorticity_cell ( nVertLevels nCells Time ) 2 o vorticity_cell state - -
-var persistent real pv_edge ( nVertLevels nEdges Time ) 2 o pv_edge state - -
var persistent real h_edge ( nVertLevels nEdges Time ) 2 o h_edge state - -
var persistent real ke ( nVertLevels nCells Time ) 2 o ke state - -
-var persistent real pv_vertex ( nVertLevels nVertices Time ) 2 o pv_vertex state - -
-var persistent real pv_cell ( nVertLevels nCells Time ) 2 o pv_cell state - -
var persistent real uReconstructX ( nVertLevels nCells Time ) 2 o uReconstructX state - -
var persistent real uReconstructY ( nVertLevels nCells Time ) 2 o uReconstructY state - -
var persistent real uReconstructZ ( nVertLevels nCells Time ) 2 o uReconstructZ state - -
var persistent real uReconstructZonal ( nVertLevels nCells Time ) 2 o uReconstructZonal state - -
var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 2 o uReconstructMeridional state - -
-
-# Other diagnostic variables: neither read nor written to any files
-var persistent real vh ( nVertLevels nEdges Time ) 2 - vh state - -
-var persistent real circulation ( nVertLevels nVertices Time ) 2 - circulation state - -
-var persistent real gradPVt ( nVertLevels nEdges Time ) 2 - gradPVt state - -
-var persistent real gradPVn ( nVertLevels nEdges Time ) 2 - gradPVn state - -
-var persistent real h_vertex ( nVertLevels nVertices Time ) 2 - h_vertex state - -
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_alpha_solver.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_alpha_solver.F 2011-04-27 16:10:24 UTC (rev 803)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_alpha_solver.F 2011-04-27 17:40:33 UTC (rev 804)
@@ -8,8 +8,8 @@
implicit none
- real (KIND=RKIND), parameter :: threshold_alpha = 1.0e-8
- integer, parameter :: nIterPerGS = 2
+ real (KIND=RKIND), parameter :: threshold_alpha = 1.0e-9
+ integer, parameter :: nIterPerGS = 2, nCheckAlpha = 25, nQuitAlpha = 1000
private
public :: init_alpha, alpha
@@ -108,7 +108,7 @@
! check the global error every so often
- if(mod(iteration,100).eq.0) then
+ if(mod(iteration,nCheckAlpha).eq.0) then
block => domain % blocklist
do while (associated(block))
call compute_alpha_error(block % state % time_levs(2) % state, block % mesh, l2_vorSmooth_error, l2_divSmooth_error, localArea)
@@ -120,12 +120,12 @@
l2_vorSmooth_error = sqrt(globalVorSmoothError/globalArea)
l2_divSmooth_error = sqrt(globalDivSmoothError/globalArea)
if(l2_vorSmooth_error.lt.threshold_alpha.and.l2_divSmooth_error.lt.threshold_alpha) converge=.true.
- if(iteration.gt.10000) converge=.true.
+ if(iteration.gt.nQuitAlpha) converge=.true.
endif
! make sure that everyone agrees on convergence
- if(mod(iteration,100).eq.0) then
+ if(mod(iteration,nCheckAlpha).eq.0) then
block => domain % blocklist
do while (associated(block))
rflag = 0.0
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_global_diagnostics.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_global_diagnostics.F 2011-04-27 16:10:24 UTC (rev 803)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_global_diagnostics.F 2011-04-27 17:40:33 UTC (rev 804)
@@ -48,7 +48,7 @@
! Step 1
! 1. Define the array to integrate, and the variable for the value to be stored in after the integration
real (kind=RKIND), dimension(:), pointer :: areaCell, dcEdge, dvEdge, areaTriangle, h_s, fCell, fEdge
- real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, pv_edge, pv_vertex, pv_cell, h_vertex, weightsOnEdge
+ real (kind=RKIND), dimension(:,:), pointer :: h, u, v, h_edge, pv, weightsOnEdge
real (kind=RKIND), dimension(:,:,:), pointer :: tracers
@@ -56,12 +56,12 @@
real (kind=RKIND), dimension(:), allocatable :: potentialEnstrophyReservior, areaEdge, h_s_edge
real (kind=RKIND), dimension(:,:), allocatable :: cellVolume, cellArea, volumeWeightedPotentialVorticity
- real (kind=RKIND), dimension(:,:), allocatable :: volumeWeightedPotentialEnstrophy, vertexVolume, volumeWeightedKineticEnergy
+ real (kind=RKIND), dimension(:,:), allocatable :: volumeWeightedPotentialEnstrophy, volumeWeightedKineticEnergy
real (kind=RKIND), dimension(:,:), allocatable :: volumeWeightedPotentialEnergy, volumeWeightedPotentialEnergyTopography
real (kind=RKIND), dimension(:,:), allocatable :: keTend_CoriolisForce, keTend_PressureGradient
real (kind=RKIND), dimension(:,:), allocatable ::peTend_DivThickness, refAreaWeightedSurfaceHeight, refAreaWeightedSurfaceHeight_edge
- real (kind=RKIND) :: sumCellVolume, sumCellArea, sumVertexVolume, sumrefAreaWeightedSurfaceHeight
+ real (kind=RKIND) :: sumCellVolume, sumCellArea, sumrefAreaWeightedSurfaceHeight
real (kind=RKIND) :: globalFluidThickness, globalPotentialVorticity, globalPotentialEnstrophy, globalEnergy
real (kind=RKIND) :: globalCoriolisEnergyTendency, globalKEPETendency, globalPotentialEnstrophyReservoir
@@ -105,10 +105,7 @@
v => state % v % array
tracers => state % tracers % array
h_edge => state % h_edge % array
- h_vertex => state % h_vertex % array
- pv_edge => state % pv_edge % array
- pv_vertex => state % pv_vertex % array
- pv_cell => state % pv_cell % array
+ pv => state % pv % array
! Step 2
! 2. Allocate the array with the correct dimensions.
@@ -116,10 +113,9 @@
allocate(cellArea(nVertLevels,nCellsSolve))
allocate(refAreaWeightedSurfaceHeight(nVertLevels,nCellsSolve))
allocate(refAreaWeightedSurfaceHeight_edge(nVertLevels,nEdgesSolve))
- allocate(volumeWeightedPotentialVorticity(nVertLevels,nVerticesSolve))
- allocate(volumeWeightedPotentialEnstrophy(nVertLevels,nVerticesSolve))
+ allocate(volumeWeightedPotentialVorticity(nVertLevels,nCellsSolve))
+ allocate(volumeWeightedPotentialEnstrophy(nVertLevels,nCellsSolve))
allocate(potentialEnstrophyReservior(nCellsSolve))
- allocate(vertexVolume(nVertLevels,nVerticesSolve))
allocate(volumeWeightedKineticEnergy(nVertLevels,nEdgesSolve))
allocate(volumeWeightedPotentialEnergy(nVertLevels,nCellsSolve))
allocate(volumeWeightedPotentialEnergyTopography(nVertLevels,nCellsSolve))
@@ -136,7 +132,6 @@
cellVolume = 0
refAreaWeightedSurfaceHeight = 0
refAreaWeightedSurfaceHeight_edge = 0
- vertexVolume = 0
cellArea = 0
averageThickness = 0
volumeWeightedPotentialVorticity = 0
@@ -159,11 +154,9 @@
cellVolume(iLevel,:) = h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)
! eg. GlobalFluidThickness bot (Sum(dA)) = Sum(cellArea)
cellArea(iLevel,:) = areaCell(1:nCellsSolve)
- volumeWeightedPotentialVorticity(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
- *h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
- volumeWeightedPotentialEnstrophy(iLevel,:) = pv_vertex(iLevel,1:nVerticesSolve) &
- *pv_vertex(iLevel,1:nVerticesSolve)*h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
- vertexVolume(iLevel,:) = h_vertex(iLevel,1:nVerticesSolve)*areaTriangle(1:nVerticesSolve)
+ volumeWeightedPotentialVorticity(iLevel,:) = pv(iLevel,1:nCellsSolve) &
+ *h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)
+ volumeWeightedPotentialEnstrophy(iLevel,:) = pv(iLevel,1:nCellsSolve)*volumeWeightedPotentialVorticity(iLevel,1:nCellsSolve)
volumeWeightedKineticEnergy(iLevel,:) = u(iLevel,1:nEdgesSolve)*u(iLevel,1:nEdgesSolve) &
*h_edge(iLevel,1:nEdgesSolve)*areaEdge(1:nEdgesSolve)*0.5
volumeWeightedPotentialEnergy(iLevel,:) = gravity*h(iLevel,1:nCellsSolve)*h(iLevel,1:nCellsSolve)*areaCell(1:nCellsSolve)*0.5
@@ -172,11 +165,11 @@
do iEdge = 1,nEdgesSolve
q = 0.0
- do j = 1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- workpv = 0.5 * (pv_edge(iLevel,iEdge) + pv_edge(iLevel,eoe))
- q = q + weightsOnEdge(j,iEdge) * u(iLevel,eoe) * workpv * h_edge(iLevel,eoe)
- end do
+ ! do j = 1,nEdgesOnEdge(iEdge)
+ ! eoe = edgesOnEdge(j,iEdge)
+ ! workpv = 0.5 * (pv_edge(iLevel,iEdge) + pv_edge(iLevel,eoe))
+ ! q = q + weightsOnEdge(j,iEdge) * u(iLevel,eoe) * workpv * h_edge(iLevel,eoe)
+ ! end do
keTend_CoriolisForce(iLevel,iEdge) = h_edge(iLevel,iEdge) * u(iLevel,iEdge) * q * areaEdge(iEdge)
iCell1 = cellsOnEdge(1,iEdge)
@@ -227,12 +220,12 @@
averageThickness(:) = (sumrefAreaWeightedSurfaceHeight/sumCellArea)-h_s(1:nCellsSolve)
! Compute Volume Weighted Averages of Potential Vorticity and Potential Enstrophy
- call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialVorticity, globalPotentialVorticity)
- call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, volumeWeightedPotentialEnstrophy, globalPotentialEnstrophy)
- call computeGlobalSum(dminfo, nVertLevels, nVerticesSolve, vertexVolume, sumVertexVolume)
+ call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialVorticity, globalPotentialVorticity)
+ call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, volumeWeightedPotentialEnstrophy, globalPotentialEnstrophy)
+ call computeGlobalSum(dminfo, nVertLevels, nCellsSolve, cellVolume, sumCellVolume)
- globalPotentialVorticity = globalPotentialVorticity/sumVertexVolume
- globalPotentialEnstrophy = globalPotentialEnstrophy/sumVertexVolume
+ globalPotentialVorticity = globalPotentialVorticity/sumCellVolume
+ globalPotentialEnstrophy = globalPotentialEnstrophy/sumCellVolume
! Compte Potential Enstrophy Reservior
potentialEnstrophyReservior(:) = areaCell(:)*fCell(:)*fCell(:)/averageThickness
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-04-27 16:10:24 UTC (rev 803)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_poisson_solver.F 2011-04-27 17:40:33 UTC (rev 804)
@@ -7,8 +7,8 @@
implicit none
- real (KIND=RKIND), parameter :: threshold = 1.0e-8
- integer, parameter :: nIterPerGS = 2
+ real (KIND=RKIND), parameter :: threshold = 1.0e-9
+ integer, parameter :: nIterPerGS = 2, nCheckPoisson = 25, nQuitPoisson = 5000
private
public :: init_poisson, poisson, derive_velocity
@@ -171,7 +171,7 @@
! check the global error every so often
- if(mod(iteration,100).eq.0) then
+ if(mod(iteration,nCheckPoisson).eq.0) then
block => domain % blocklist
do while (associated(block))
call compute_poisson_error(block % state % time_levs(2) % state, block % mesh, l2_psi_error, l2_chi_error, localArea)
@@ -183,12 +183,12 @@
l2_psi_error = sqrt(globalPsiError/globalArea)
l2_chi_error = sqrt(globalChiError/globalArea)
if(l2_psi_error.lt.threshold.and.l2_chi_error.lt.threshold) converge=.true.
- if(iteration.gt.10000) converge=.true.
+ if(iteration.gt.nQuitPoisson) converge=.true.
endif
! make sure that everyone agrees on convergence
- if(mod(iteration,100).eq.0) then
+ if(mod(iteration,nCheckPoisson).eq.0) then
block => domain % blocklist
do while (associated(block))
rflag = 0.0
Modified: branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F
===================================================================
--- branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-04-27 16:10:24 UTC (rev 803)
+++ branches/pv_based_swm/mpas/src/core_pvsw/module_time_integration.F 2011-04-27 17:40:33 UTC (rev 804)
@@ -7,6 +7,7 @@
use poisson_solver
use alpha_solver
use dmpar
+ use timer
contains
@@ -70,6 +71,7 @@
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
+ call timer_start("rk4 init")
block => domain % blocklist
call allocate_state(provis, &
block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
@@ -116,7 +118,9 @@
rk_substep_weights(3) = dt
rk_substep_weights(4) = 0.
+ call timer_stop("rk4 init")
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN RK loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -126,10 +130,7 @@
! block => domain % blocklist
! do while (associated(block))
- ! call dmpar_exch_halo_field2dReal(domain % dminfo, provis % pv_edge % array(:,:), &
- ! block % mesh % nVertLevels, block % mesh % nEdges, &
- ! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
-
+ !
! if (config_h_mom_eddy_visc4 > 0.0) then
! call dmpar_exch_halo_field2dReal(domain % dminfo, block % state % time_levs(2) % state % divergence % array(:,:), &
! block % mesh % nVertLevels, block % mesh % nCells, &
@@ -155,11 +156,14 @@
! --- compute smoothed vorticity for alpha closure
if(config_alpha_closure) then
+ call timer_start("alpha closure")
call alpha(domain)
+ call timer_stop("alpha closure")
endif
! --- invert elliptic equations for psi and chi and then derive normal and tangent velocity
+ call timer_start("poisson")
call poisson(domain)
call derive_velocity(domain)
block => domain % blocklist
@@ -168,9 +172,11 @@
provis % v % array(:,:) = block % state % time_levs(2) % state % v % array(:,:)
block => block % next
end do
+ call timer_stop("poisson")
! --- compute tendencies
+ call timer_start("tendencies")
block => domain % blocklist
do while (associated(block))
call compute_tend(block % tend, provis, block % mesh)
@@ -178,14 +184,13 @@
call enforce_boundaryEdge(block % tend, block % mesh)
block => block % next
end do
+ call timer_stop("tendencies")
! --- update halos for prognostic variables
+ call timer_start("prognostic halo")
block => domain % blocklist
do while (associated(block))
- ! call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % u % array(:,:), &
- ! block % mesh % nVertLevels, block % mesh % nEdges, &
- ! block % parinfo % edgesToSend, block % parinfo % edgesToRecv)
call dmpar_exch_halo_field2dReal(domain % dminfo, block % tend % h % array(:,:), &
block % mesh % nVertLevels, block % mesh % nCells, &
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
@@ -194,9 +199,11 @@
block % parinfo % cellsToSend, block % parinfo % cellsToRecv)
block => block % next
end do
+ call timer_stop("prognostic halo")
! --- compute next substep state
+ call timer_start("update provisional")
if (rk_step < 4) then
block => domain % blocklist
do while (associated(block))
@@ -216,9 +223,11 @@
block => block % next
end do
end if
+ call timer_stop("update provisional")
!--- accumulate update (for RK4)
+ call timer_start("increment prognostic")
block => domain % blocklist
do while (associated(block))
block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(2) % state % h % array(:,:) &
@@ -232,6 +241,7 @@
end do
block => block % next
end do
+ call timer_stop("increment prognostic")
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -243,6 +253,7 @@
! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
!
+ call timer_start("rk4 cleanup")
block => domain % blocklist
do while (associated(block))
do iCell=1,block % mesh % nCells
@@ -286,6 +297,7 @@
block => block % next
end do
+ call timer_stop("rk4 cleanup")
call deallocate_state(provis)
@@ -309,19 +321,18 @@
type (mesh_type), intent(in) :: grid
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
- real (kind=RKIND) :: flux, vorticity_abs, workpv, q, upstream_bias
+ real (kind=RKIND) :: flux, q, upstream_bias
integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
- real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- circulation, vorticity, ke, pv_edge, divergence, h_vertex
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnVertex, verticesOnEdge, edgesOnCell, edgesOnEdge, edgesOnVertex
+ real (kind=RKIND), dimension(:), pointer :: h_s, dvEdge, dcEdge, areaCell
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, h_edge, h, u, v, tend_h, tend_u, &
+ ke, divergence
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, edgesOnEdge
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
real (kind=RKIND) :: r, u_diffusion
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_divergence
real (kind=RKIND), allocatable, dimension(:,:) :: delsq_u
- real (kind=RKIND), allocatable, dimension(:,:) :: delsq_circulation, delsq_vorticity
real (kind=RKIND), dimension(:,:), pointer :: u_src
real (kind=RKIND), parameter :: rho_ref = 1000.0
@@ -331,30 +342,19 @@
u => s % u % array
v => s % v % array
h_edge => s % h_edge % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
divergence => s % divergence % array
ke => s % ke % array
- pv_edge => s % pv_edge % array
- vh => s % vh % array
weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
nEdgesOnCell => grid % nEdgesOnCell % array
edgesOnCell => grid % edgesOnCell % array
nEdgesOnEdge => grid % nEdgesOnEdge % array
edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
dcEdge => grid % dcEdge % array
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
h_s => grid % h_s % array
- fVertex => grid % fVertex % array
- fEdge => grid % fEdge % array
tend_h => tend % h % array
tend_u => tend % u % array
@@ -687,16 +687,15 @@
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
- real (kind=RKIND) :: flux, vorticity_abs, workpv
+ real (kind=RKIND) :: flux
integer :: nCells, nEdges, nVertices, nVertLevels
- real (kind=RKIND), dimension(:), pointer :: h_s, fVertex, fEdge, dvEdge, dcEdge, areaCell, areaTriangle
- real (kind=RKIND), dimension(:,:), pointer :: vh, weightsOnEdge, kiteAreasOnVertex, h_edge, h, u, v, tend_h, tend_u, &
- 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
+ real (kind=RKIND), dimension(:), pointer :: h_s, dvEdge, dcEdge, areaCell
+ real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge, h_edge, h, u, v, tend_h, tend_u, &
+ ke, divergence
+ integer, dimension(:,:), pointer :: cellsOnEdge, edgesOnCell, edgesOnEdge, boundaryEdge, boundaryCell
integer, dimension(:), pointer :: nEdgesOnCell, nEdgesOnEdge
- real (kind=RKIND) :: r, h1, h2
+ real (kind=RKIND) :: r
real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
real (kind=RKIND) :: coef_3rd_order
@@ -704,39 +703,22 @@
h => s % h % array
u => s % u % array
v => s % v % array
- vh => s % vh % array
h_edge => s % h_edge % array
- h_vertex => s % h_vertex % array
tend_h => s % h % array
tend_u => s % u % array
- circulation => s % circulation % array
- vorticity => s % vorticity % array
divergence => s % divergence % array
ke => s % ke % array
- pv_edge => s % pv_edge % array
- pv_vertex => s % pv_vertex % array
- pv_cell => s % pv_cell % array
- vorticity_cell => s % vorticity_cell % array
- gradPVn => s % gradPVn % array
- gradPVt => s % gradPVt % array
weightsOnEdge => grid % weightsOnEdge % array
- kiteAreasOnVertex => grid % kiteAreasOnVertex % array
cellsOnEdge => grid % cellsOnEdge % array
- cellsOnVertex => grid % cellsOnVertex % array
- verticesOnEdge => grid % verticesOnEdge % array
nEdgesOnCell => grid % nEdgesOnCell % array
edgesOnCell => grid % edgesOnCell % array
nEdgesOnEdge => grid % nEdgesOnEdge % array
edgesOnEdge => grid % edgesOnEdge % array
- edgesOnVertex => grid % edgesOnVertex % array
dcEdge => grid % dcEdge % array
dvEdge => grid % dvEdge % array
areaCell => grid % areaCell % array
- areaTriangle => grid % areaTriangle % array
h_s => grid % h_s % array
- fVertex => grid % fVertex % array
- fEdge => grid % fEdge % array
deriv_two => grid % deriv_two % array
nCells = grid % nCells
@@ -763,162 +745,24 @@
enddo
!
- ! Compute height on cell edges at velocity locations
- ! Namelist options control the order of accuracy of the reconstructed h_edge value
- !
-
- coef_3rd_order = 0.
- if (config_thickness_adv_order == 3) coef_3rd_order = 1.0
- if (config_thickness_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
-
- if (config_thickness_adv_order == 2) then
-
- 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
- h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
- end do
- end if
- end do
-
- else if (config_thickness_adv_order == 3) then
-
- 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
-
- !-- 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) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(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) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- !-- if u > 0:
- if (u(k,iEdge) > 0) then
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(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
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
- +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12.
- end if
-
- end do ! do k
- end if ! if (cell1 <=
- end do ! do iEdge
-
- else if (config_thickness_adv_order == 4) then
-
- 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
-
- !-- 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) * h(k,cell1)
- d2fdx2_cell2 = deriv_two(1,2,iEdge) * h(k,cell2)
-
- !-- all edges of cell 1
- do i=1, grid % nEdgesOnCell % array (cell1)
- d2fdx2_cell1 = d2fdx2_cell1 + &
- deriv_two(i+1,1,iEdge) * h(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) * h(k,grid % CellsOnCell % array (i,cell2))
- end do
-
- endif
-
- h_edge(k,iEdge) = &
- 0.5*(h(k,cell1) + h(k,cell2)) &
- -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
-
- end do ! do k
- end if ! if (cell1 <=
- end do ! do iEdge
-
- endif ! if(config_thickness_adv_order == 2)
-
- !
! set the velocity in the nEdges+1 slot to zero, this is a dummy address
! used to when reading for edges that do not exist
!
u(:,nEdges+1) = 0.0
!
- ! Compute circulation and relative vorticity at each vertex
- !
- circulation(:,:) = 0.0
- do iEdge=1,nEdges
- do k=1,nVertLevels
- circulation(k,verticesOnEdge(1,iEdge)) = circulation(k,verticesOnEdge(1,iEdge)) - dcEdge(iEdge) * u(k,iEdge)
- circulation(k,verticesOnEdge(2,iEdge)) = circulation(k,verticesOnEdge(2,iEdge)) + dcEdge(iEdge) * u(k,iEdge)
- end do
- end do
- do iVertex=1,nVertices
- do k=1,nVertLevels
- vorticity(k,iVertex) = circulation(k,iVertex) / areaTriangle(iVertex)
- end do
- end do
-
-
- !
! Compute the divergence at each cell center
!
divergence(:,:) = 0.0
+ h_edge(:,:) = 0.0
do iEdge=1,nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= nCells) then
- do k=1,nVertLevels
- divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
- enddo
- endif
- if(cell2 <= nCells) then
- do k=1,nVertLevels
- divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
- enddo
- end if
+ do k=1,nVertLevels
+ divergence(k,cell1) = divergence(k,cell1) + u(k,iEdge)*dvEdge(iEdge)
+ divergence(k,cell2) = divergence(k,cell2) - u(k,iEdge)*dvEdge(iEdge)
+ h_edge(k,iEdge) = 0.5 * (h(k,cell1) + h(k,cell2))
+ enddo
end do
do iCell = 1,nCells
r = 1.0 / areaCell(iCell)
@@ -956,141 +800,6 @@
end do
end do
-#ifdef NCAR_FORMULATION
- !
- ! Compute mass fluxes tangential to each edge (i.e., through the faces of dual grid cells)
- !
- vh(:,:) = 0.0
- do iEdge=1,grid % nEdgesSolve
- do j=1,nEdgesOnEdge(iEdge)
- eoe = edgesOnEdge(j,iEdge)
- do k=1,nVertLevels
- vh(k,iEdge) = vh(k,iEdge) + weightsOnEdge(j,iEdge) * u(k,eoe) * h_edge(k,eoe)
- end do
- end do
- 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 )
- !
- do iVertex = 1,nVertices
- do k=1,nVertLevels
- h_vertex(k,iVertex) = 0.0
- do i=1,grid % vertexDegree
- h_vertex(k,iVertex) = h_vertex(k,iVertex) + h(k,cellsOnVertex(i,iVertex)) * kiteAreasOnVertex(i,iVertex)
- end do
- h_vertex(k,iVertex) = h_vertex(k,iVertex) / areaTriangle(iVertex)
-
- pv_vertex(k,iVertex) = (fVertex(iVertex) + vorticity(k,iVertex)) / h_vertex(k,iVertex)
- end do
- end do
-
-
- !
- ! Compute gradient of PV in the tangent direction
- ! ( this computes gradPVt at all edges bounding real cells and distance-1 ghost cells )
- !
- do iEdge = 1,nEdges
- do k = 1,nVertLevels
- gradPVt(k,iEdge) = (pv_vertex(k,verticesOnEdge(2,iEdge)) - pv_vertex(k,verticesOnEdge(1,iEdge))) / &
- dvEdge(iEdge)
- enddo
- enddo
-
- !
- ! Compute pv at the edges
- ! ( this computes pv_edge at all edges bounding real cells )
- !
- pv_edge(:,:) = 0.0
- 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
- end do
- end do
-
-
- !
- ! Modify PV edge with upstream bias.
- !
- 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)
- enddo
- enddo
-
-
- !
- ! Compute pv at cell centers
- ! ( this computes pv_cell for all real cells and distance-1 ghost cells )
- !
- pv_cell(:,:) = 0.0
- vorticity_cell(:,:) = 0.0
- 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
-
-
- !
- ! Compute gradient of PV in normal direction
- ! ( this computes gradPVn for all edges bounding real cells )
- !
- gradPVn(:,:) = 0.0
- do iEdge = 1,nEdges
- if( cellsOnEdge(1,iEdge) <= nCells .and. cellsOnEdge(2,iEdge) <= nCells) then
- do k = 1,nVertLevels
- gradPVn(k,iEdge) = (pv_cell(k,cellsOnEdge(2,iEdge)) - pv_cell(k,cellsOnEdge(1,iEdge))) / &
- dcEdge(iEdge)
- enddo
- endif
- enddo
-
- ! Modify PV edge with upstream bias.
- !
- 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)
- enddo
- enddo
-
- !
- ! set pv_edge = fEdge / h_edge at boundary points
- !
- ! if (maxval(boundaryEdge).ge.0) then
- ! do iEdge = 1,nEdges
- ! cell1 = cellsOnEdge(1,iEdge)
- ! cell2 = cellsOnEdge(2,iEdge)
- ! do k = 1,nVertLevels
- ! if(boundaryEdge(k,iEdge).eq.1) then
- ! v(k,iEdge) = 0.0
- ! if(cell1.gt.0) then
- ! h1 = h(k,cell1)
- ! pv_edge(k,iEdge) = fEdge(iEdge) / h1
- ! h_edge(k,iEdge) = h1
- ! else
- ! h2 = h(k,cell2)
- ! pv_edge(k,iEdge) = fEdge(iEdge) / h2
- ! h_edge(k,iEdge) = h2
- ! endif
- ! endif
- ! enddo
- ! enddo
- ! endif
-
-
end subroutine compute_solve_diagnostics
</font>
</pre>