<p><b>dwj07@fsu.edu</b> 2012-10-22 15:24:03 -0600 (Mon, 22 Oct 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Adding another optional input argument.<br>
This argument represents the missing factor between ocean and atmosphere for vertical divergences.<br>
<br>
If it's not present on input it defaults to 1.0.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_mono.F 2012-10-22 21:21:30 UTC (rev 2247)
+++ branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_mono.F 2012-10-22 21:24:03 UTC (rev 2248)
@@ -44,7 +44,7 @@
!> Both horizontal and vertical.
!
!-----------------------------------------------------------------------
- subroutine mpas_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend, maxLevelCell_in, maxLevelEdgeTop_in, highOrderAdvectionMask_in)!{{{
+ subroutine mpas_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend, maxLevelCell_in, maxLevelEdgeTop_in, highOrderAdvectionMask_in, verticalDivergenceFactor_in)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -63,7 +63,8 @@
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend !< Input/Output: Tracer tendency
integer, dimension(:), pointer, optional :: maxLevelCell_in !< Input - optional: Index to max level at cell center
integer, dimension(:), pointer, optional :: maxLevelEdgeTop_in !< Input - optional: Index to max level at edge with non-land cells on both sides
- integer, dimension(:,:), pointer, optional :: highOrderAdvectionMask_in !< Input- optional: Mask for high order advection
+ integer, dimension(:,:), pointer, optional :: highOrderAdvectionMask_in !< Input - optional: Mask for high order advection
+ real (kind=RKIND), dimension(:), pointer, optional :: verticalDivergenceFactor_in !< Input - optional: Denominator for vertical divergence. Given as an inverse which can be multiplied.
integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
integer :: nVertLevels, num_tracers, nCells, nEdges, nCellsSolve
@@ -74,12 +75,12 @@
real (kind=RKIND) :: flux, tracer_weight, invDvEdge, invAreaCell1, invAreaCell2
real (kind=RKIND) :: cur_max, cur_min, new_max, new_min
real (kind=RKIND) :: verticalWeightK, verticalWeightKm1
- real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, verticalDivergenceFactor
real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_2nd, adv_coefs_3rd
real (kind=RKIND), dimension(:,:), pointer :: tracer_cur, tracer_new, upwind_tendency, inv_h_new, tracer_max, tracer_min
real (kind=RKIND), dimension(:,:), pointer :: flux_incoming, flux_outgoing, high_order_horiz_flux, high_order_vert_flux
- real (kind=RKIND), parameter :: eps = 1.e-10
+ real (kind=RKIND), parameter :: eps = 1.e-10_RKIND
type (field2dReal), pointer :: high_order_horiz_flux_field
@@ -126,7 +127,13 @@
else
highOrderAdvectionMask = 1
end if
+ end if
+ if(present(verticalDivergenceFactor_in)) then
+ verticalDivergenceFactor => verticalDivergenceFactor_in
+ else
+ allocate(verticalDivergenceFactor(nVertLevels))
+ verticalDivergenceFactor = 1.0_RKIND
end if
! Setup high order horizontal flux field
@@ -156,7 +163,7 @@
do iCell = 1, nCells
do k=1, maxLevelCell(iCell)
- inv_h_new(k, iCell) = 1.0 / (h(k, iCell) + dt * tend_h(k, iCell))
+ inv_h_new(k, iCell) = 1.0_RKIND / (h(k, iCell) + dt * tend_h(k, iCell))
end do
end do
@@ -166,17 +173,17 @@
do iCell = 1, nCells
do k=1, maxLevelCell(iCell)
tracer_cur(k,iCell) = tracers(iTracer,k,iCell)
- upwind_tendency(k, iCell) = 0.0
+ upwind_tendency(k, iCell) = 0.0_RKIND
!tracer_new is supposed to be the "new" tracer state. This allows bounds checks.
if (config_check_monotonicity) then
- tracer_new(k,iCell) = 0.0
+ tracer_new(k,iCell) = 0.0_RKIND
end if
end do ! k loop
end do ! iCell loop
- high_order_vert_flux = 0.0
- high_order_horiz_flux = 0.0
+ high_order_vert_flux = 0.0_RKIND
+ high_order_horiz_flux = 0.0_RKIND
! Compute the high order vertical flux. Also determine bounds on tracer_cur.
do iCell = 1, nCells
@@ -231,7 +238,7 @@
! Compute 2nd order fluxes where needed.
do k = 1, maxLevelEdgeTop(iEdge)
- tracer_weight = iand(highOrderAdvectionMask(k, iEdge)+1, 1) * (dvEdge(iEdge) * 0.5) * uh(k, iEdge)
+ tracer_weight = iand(highOrderAdvectionMask(k, iEdge)+1, 1) * (dvEdge(iEdge) * 0.5_RKIND) * uh(k, iEdge)
high_order_horiz_flux(k, iEdge) = high_order_horiz_flux(k, iedge) + tracer_weight * (tracer_cur(k, cell1) + tracer_cur(k, cell2))
end do ! k loop
@@ -240,7 +247,7 @@
do i = 1, nAdvCellsForEdge(iEdge)
iCell = advCellsForEdge(i,iEdge)
do k = 1, maxLevelCell(iCell)
- tracer_weight = highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
+ tracer_weight = highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
tracer_weight = uh(k,iEdge)*tracer_weight
high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) + tracer_weight * tracer_cur(k,iCell)
@@ -258,9 +265,9 @@
do k = 2, maxLevelCell(iCell)
! dwj 02/03/12 and Atmosphere are different in vertical
if(config_dzdk_positive) then
- flux_upwind = max(0.,w(k,iCell))*tracer_cur(k-1,iCell) + min(0.,w(k,iCell))*tracer_cur(k,iCell)
+ flux_upwind = max(0.0_RKIND,w(k,iCell))*tracer_cur(k-1,iCell) + min(0.0_RKIND,w(k,iCell))*tracer_cur(k,iCell)
else
- flux_upwind = min(0.,w(k,iCell))*tracer_cur(k-1,iCell) + max(0.,w(k,iCell))*tracer_cur(k,iCell)
+ flux_upwind = min(0.0_RKIND,w(k,iCell))*tracer_cur(k-1,iCell) + max(0.0_RKIND,w(k,iCell))*tracer_cur(k,iCell)
end if
upwind_tendency(k-1,iCell) = upwind_tendency(k-1,iCell) + flux_upwind
upwind_tendency(k ,iCell) = upwind_tendency(k ,iCell) - flux_upwind
@@ -274,11 +281,11 @@
do k = 1, maxLevelCell(iCell)
! dwj 02/03/12 and Atmosphere are different in vertical
if(config_dzdk_positive) then
- flux_incoming (k,iCell) = -(min(0.,high_order_vert_flux(k+1,iCell))-max(0.,high_order_vert_flux(k,iCell)))
- flux_outgoing(k,iCell) = -(max(0.,high_order_vert_flux(k+1,iCell))-min(0.,high_order_vert_flux(k,iCell)))
+ flux_incoming (k,iCell) = -(min(0.0_RKIND,high_order_vert_flux(k+1,iCell))-max(0.0_RKIND,high_order_vert_flux(k,iCell)))
+ flux_outgoing(k,iCell) = -(max(0.0_RKIND,high_order_vert_flux(k+1,iCell))-min(0.0_RKIND,high_order_vert_flux(k,iCell)))
else
- flux_incoming (k, iCell) = max(0.0, high_order_vert_flux(k+1, iCell)) - min(0.0, high_order_vert_flux(k, iCell))
- flux_outgoing(k, iCell) = min(0.0, high_order_vert_flux(k+1, iCell)) - max(0.0, high_order_vert_flux(k, iCell))
+ flux_incoming (k, iCell) = max(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - min(0.0_RKIND, high_order_vert_flux(k, iCell))
+ flux_outgoing(k, iCell) = min(0.0_RKIND, high_order_vert_flux(k+1, iCell)) - max(0.0_RKIND, high_order_vert_flux(k, iCell))
end if
end do ! k Loop
end do ! iCell Loop
@@ -291,8 +298,8 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
+ invAreaCell1 = 1.0_RKIND / areaCell(cell1)
+ invAreaCell2 = 1.0_RKIND / areaCell(cell2)
do k = 1, maxLevelEdgeTop(iEdge)
flux_upwind = dvEdge(iEdge) * (max(0.,uh(k,iEdge))*tracer_cur(k,cell1) + min(0.,uh(k,iEdge))*tracer_cur(k,cell2))
@@ -319,10 +326,10 @@
tracer_upwind_new = (tracer_cur(k,iCell)*h(k,iCell) + dt*upwind_tendency(k,iCell)) * inv_h_new(k,iCell)
scale_factor = (tracer_max(k,iCell)-tracer_upwind_new)/(tracer_max_new-tracer_upwind_new+eps)
- flux_incoming(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
+ flux_incoming(k,iCell) = min( 1.0_RKIND, max( 0.0_RKIND, scale_factor) )
scale_factor = (tracer_upwind_new-tracer_min(k,iCell))/(tracer_upwind_new-tracer_min_new+eps)
- flux_outgoing(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
+ flux_outgoing(k,iCell) = min( 1.0_RKIND, max( 0.0_RKIND, scale_factor) )
end do ! k loop
end do ! iCell loop
@@ -359,8 +366,8 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
+ invAreaCell1 = 1.0_RKIND / areaCell(cell1)
+ invAreaCell2 = 1.0_RKIND / areaCell(cell2)
do k=1,maxLevelEdgeTop(iEdge)
tend(iTracer, k, cell1) = tend(iTracer, k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
tend(iTracer, k, cell2) = tend(iTracer, k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
@@ -376,11 +383,11 @@
! Accumulate the scaled high order vertical tendencies, and the upwind tendencies
do iCell = 1, nCellsSolve
do k = 1,maxLevelCell(iCell)
- tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + verticalDivergenceFactor(k) * (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
if (config_check_monotonicity) then
!tracer_new holds a tendency for now. Only for a check on monotonicity
- tracer_new(k, iCell) = tracer_new(k, iCell) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
+ tracer_new(k, iCell) = tracer_new(k, iCell) + verticalDivergenceFactor(k) * (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
!tracer_new is now the new state of the tracer. Only for a check on monotonicity
tracer_new(k, iCell) = (tracer_cur(k, iCell)*h(k, iCell) + dt * tracer_new(k, iCell)) * inv_h_new(k, iCell)
@@ -427,6 +434,10 @@
if(.not.present(highOrderAdvectionMask_in)) then
deallocate(highOrderAdvectionMask)
end if
+
+ if(.not.present(verticalDivergenceFactor_in)) then
+ deallocate(verticalDivergenceFactor)
+ end if
end subroutine mpas_tracer_advection_mono_tend!}}}
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
@@ -455,7 +466,7 @@
if ( config_horiz_tracer_adv_order == 3) then
coef_3rd_order = config_coef_3rd_order
else if(config_horiz_tracer_adv_order == 2 .or. config_horiz_tracer_adv_order == 4) then
- coef_3rd_order = 0.0
+ coef_3rd_order = 0.0_RKIND
end if
if (config_vert_tracer_adv_order == 3) then
Modified: branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_std.F
===================================================================
--- branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_std.F 2012-10-22 21:21:30 UTC (rev 2247)
+++ branches/ocean_projects/shared_advection/src/operators/mpas_tracer_advection_std.F 2012-10-22 21:24:03 UTC (rev 2248)
@@ -46,7 +46,7 @@
!> Both horizontal and vertical.
!
!-----------------------------------------------------------------------
- subroutine mpas_tracer_advection_std_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend, maxLevelCell_in, maxLevelEdgeTop_in, highOrderAdvectionMask_in)!{{{
+ subroutine mpas_tracer_advection_std_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend, maxLevelCell_in, maxLevelEdgeTop_in, highOrderAdvectionMask_in, verticalDivergenceFactor_in)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -66,6 +66,7 @@
integer, dimension(:), pointer, optional :: maxLevelCell_in !< Input - optional: Index to max level at cell center
integer, dimension(:), pointer, optional :: maxLevelEdgeTop_in !< Input - optional: Index to max level at edge with non-land cells on both sides
integer, dimension(:,:), pointer, optional :: highOrderAdvectionMask_in !< Input- optional: Mask for high order advection
+ real (kind=RKIND), dimension(:), pointer, optional :: verticalDivergenceFactor_in !< Input - optional: Denominator for vertical divergence. Given as an inverse which can be multiplied.
integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
integer :: nVertLevels, num_tracers, nCells, nEdges, nCellsSolve
@@ -74,11 +75,11 @@
real (kind=RKIND) :: flux, tracer_weight, invDvEdge, invAreaCell1, invAreaCell2
real (kind=RKIND) :: verticalWeightK, verticalWeightKm1
- real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell, verticalDivergenceFactor
real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_2nd, adv_coefs_3rd
real (kind=RKIND), dimension(:,:), pointer :: high_order_horiz_flux, high_order_vert_flux, tracer_cur
- real (kind=RKIND), parameter :: eps = 1.e-10
+ real (kind=RKIND), parameter :: eps = 1.e-10_RKIND
type (field2dReal), pointer :: high_order_horiz_flux_field
@@ -125,7 +126,13 @@
else
highOrderAdvectionMask = 1
end if
+ end if
+ if(present(verticalDivergenceFactor_in)) then
+ verticalDivergenceFactor => verticalDivergenceFactor_in
+ else
+ allocate(verticalDivergenceFactor(nVertLevels))
+ verticalDivergenceFactor = 1.0_RKIND
end if
! Setup high order horizontal flux field
@@ -149,8 +156,8 @@
! Loop over tracers. One tracer is advected at a time. It is copied into a temporary array in order to improve locality
do iTracer = 1, num_tracers
! Initialize variables for use in this iTracer iteration
- high_order_vert_flux = 0.0
- high_order_horiz_flux = 0.0
+ high_order_vert_flux = 0.0_RKIND
+ high_order_horiz_flux = 0.0_RKIND
tracer_cur(:,:) = tracers(iTracer, :, :)
! Compute the high order vertical flux. Also determine bounds on tracer_cur.
@@ -187,7 +194,7 @@
! Compute 2nd order fluxes where needed.
do k = 1, maxLevelEdgeTop(iEdge)
- tracer_weight = iand(highOrderAdvectionMask(k, iEdge)+1, 1) * (dvEdge(iEdge) * 0.5) * uh(k, iEdge)
+ tracer_weight = iand(highOrderAdvectionMask(k, iEdge)+1, 1) * (dvEdge(iEdge) * 0.5_RKIND) * uh(k, iEdge)
high_order_horiz_flux(k, iEdge) = high_order_horiz_flux(k, iedge) + tracer_weight * (tracer_cur(k, cell1) + tracer_cur(k, cell2))
end do ! k loop
@@ -196,7 +203,7 @@
do i = 1, nAdvCellsForEdge(iEdge) * highOrderHorizOn
iCell = advCellsForEdge(i,iEdge)
do k = 1, maxLevelCell(iCell)
- tracer_weight = highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
+ tracer_weight = highOrderAdvectionMask(k, iEdge) * (adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0_RKIND,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
tracer_weight = uh(k,iEdge)*tracer_weight
high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) + tracer_weight * tracer_cur(k,iCell)
@@ -209,18 +216,18 @@
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
+ invAreaCell1 = 1.0_RKIND / areaCell(cell1)
+ invAreaCell2 = 1.0_RKIND / areaCell(cell2)
do k=1,maxLevelEdgeTop(iEdge)
tend(iTracer, k, cell1) = tend(iTracer, k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
tend(iTracer, k, cell2) = tend(iTracer, k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
end do ! k loop
end do ! iEdge loop
- ! Accumulate the scaled high order vertical tendencies, and the upwind tendencies
+ ! Accumulate the scaled high order vertical tendencies.
do iCell = 1, nCellsSolve
do k = 1,maxLevelCell(iCell)
- tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell))
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + verticalDivergenceFactor(k) * (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell))
end do ! k loop
end do ! iCell loop
@@ -242,6 +249,10 @@
if(.not.present(highOrderAdvectionMask_in)) then
deallocate(highOrderAdvectionMask)
end if
+
+ if(.not.present(verticalDivergenceFactor_in)) then
+ deallocate(verticalDivergenceFactor)
+ end if
end subroutine mpas_tracer_advection_std_tend!}}}
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
</font>
</pre>