<p><b>dwj07@fsu.edu</b> 2012-03-08 20:20:50 -0700 (Thu, 08 Mar 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Setting up weights for vertical advection to support other vertical grids.<br>
Cleaning up interfaces of new shared advection routines.<br>
</p><hr noshade><pre><font color="gray">
Property changes on: branches/ocean_projects/monotonic_advection
___________________________________________________________________
Added: svn:mergeinfo
+ /trunk/mpas:1508
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/Registry 2012-03-08 23:21:56 UTC (rev 1609)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/Registry 2012-03-09 03:20:50 UTC (rev 1610)
@@ -178,6 +178,7 @@
var persistent integer maxLevelVertexTop ( nVertices ) 0 - maxLevelVertexTop mesh - -
var persistent integer maxLevelVertexBot ( nVertices ) 0 - maxLevelVertexBot mesh - -
var persistent real referenceBottomDepth ( nVertLevels ) 0 iro referenceBottomDepth mesh - -
+var persistent real referenceBottomDepthTopOfCell ( nVertLevelsP1 ) 0 - referenceBottomDepthTopOfCell mesh - -
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
var persistent real hMeanTopZLevel ( nVertLevels ) 0 - hMeanTopZLevel mesh - -
var persistent real hRatioZLevelK ( nVertLevels ) 0 - hRatioZLevelK mesh - -
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_mpas_core.F 2012-03-08 23:21:56 UTC (rev 1609)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_mpas_core.F 2012-03-09 03:20:50 UTC (rev 1610)
@@ -514,7 +514,8 @@
integer :: iTracer, cell, cell1, cell2
real (kind=RKIND) :: uhSum, hSum, hEdge1
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth, hMeanTopZLevel, hRatioZLevelK, hRatioZLevelKm1
+ real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth, referenceBottomDepthTopOfCell, hMeanTopZLevel, hRatioZLevelK, hRatioZLevelKm1
+
real (kind=RKIND), dimension(:,:), pointer :: h
integer :: nVertLevels
@@ -524,6 +525,7 @@
h => block % state % time_levs(1) % state % h % array
referenceBottomDepth => block % mesh % referenceBottomDepth % array
+ referenceBottomDepthTopOfCell => block % mesh % referenceBottomDepthTopOfCell % array
nVertLevels = block % mesh % nVertLevels
hMeanTopZLevel => block % mesh % hMeanTopZLevel % array
hRatioZLevelK => block % mesh % hRatioZLevelK % array
@@ -545,6 +547,11 @@
hRatioZLevelKm1(k) = 0.5*block % mesh % hZLevel % array(k-1)/hMeanTopZLevel(k)
end do
+ ! TopOfCell needed where zero depth for the very top may be referenced.
+ referenceBottomDepthTopOfCell(1) = 0.0
+ do k = 1,nVertLevels
+ referenceBottomDepthTopOfCell(k+1) = referenceBottomDepth(k)
+ end do
! Compute barotropic velocity at first timestep
! This is only done upon start-up.
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tendency.F 2012-03-08 23:21:56 UTC (rev 1609)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tendency.F 2012-03-09 03:20:50 UTC (rev 1610)
@@ -353,13 +353,9 @@
! and then change maxLevelEdgeTop to maxLevelEdgeBot in the following section.
! tracer_edge at the boundary will also need to be defined for flux boundaries.
- call mpas_timer_start("adv", .false., tracerHadvTimer)
! Monotonoic Advection, or standard advection
- call mpas_tracer_advection_tend(tracers, uh, wTop, h, dt, hMeanTopZLevel, hRatioZLevelKm1, hRatioZLevelK, grid, tend_h, tend_tr)
-
- ! Only standadrd advection
-! call ocn_tracer_hadv_tend(grid, u, h_edge, tracers, tend_tr, err)
-! call ocn_tracer_vadv_tend(grid, h, wtop, tracers, tend_tr, err)
+ call mpas_timer_start("adv", .false., tracerHadvTimer)
+ call mpas_tracer_advection_tend(tracers, uh, wTop, h, h, dt, grid, tend_h, tend_tr)
call mpas_timer_stop("adv", tracerHadvTimer)
!
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_split.F 2012-03-08 23:21:56 UTC (rev 1609)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_split.F 2012-03-09 03:20:50 UTC (rev 1610)
@@ -432,7 +432,7 @@
sshEdge = 0.5 * (block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell1) &
+ block % state % time_levs(oldBtrSubcycleTime) % state % sshSubcycle % array(cell2) )
- hSum = sshEdge + block % mesh % referenceBottomDepth % array (max(block % mesh % maxLevelEdgeTop % array(iEdge),1))
+ hSum = sshEdge + block % mesh % referenceBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
flux = ((1.0-config_btr_gam1_uWt1) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ config_btr_gam1_uWt1 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
@@ -546,7 +546,7 @@
+ config_btr_gam2_SSHWt1 *block % state % time_levs(newBtrSubcycleTime) % state % sshSubcycle % array(cell2)
sshEdge = 0.5 * (sshCell1 + sshCell2)
- hSum = sshEdge + block % mesh % referenceBottomDepth % array (block % mesh % maxLevelEdgeTop % array(iEdge))
+ hSum = sshEdge + block % mesh % referenceBottomDepthTopOfCell % array (block % mesh % maxLevelEdgeTop % array(iEdge)+1)
flux = ((1.0-config_btr_gam3_uWt2) * block % state % time_levs(oldBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge) &
+ config_btr_gam3_uWt2 * block % state % time_levs(newBtrSubcycleTime) % state % uBtrSubcycle % array(iEdge)) &
Modified: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection.F 2012-03-08 23:21:56 UTC (rev 1609)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection.F 2012-03-09 03:20:50 UTC (rev 1610)
@@ -176,26 +176,21 @@
end subroutine mpas_tracer_advection_coefficients!}}}
-! subroutine mpas_tracer_advection_tend(tracers, uh, w, zDistance, zWeightK, zWeightKm1, grid, dt, dminfo, cellsToSend, cellsToRecv, tend)
- subroutine mpas_tracer_advection_tend(tracers, uh, w, h, dt, zDistance, zWeightK, zWeightKm1, grid, tend_h, tend)!{{{
+ subroutine mpas_tracer_advection_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend)!{{{
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
real (kind=RKIND), dimension(:,:), intent(in) :: uh, w, h
- real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !< Input: Distance between vertical interfaces of a cell
real (kind=RKIND), intent(in) :: dt
type (mesh_type), intent(in) :: grid
real (kind=RKIND), dimension(:,:), intent(in) :: tend_h
-! real (kind=RKIND) :: dt
-! type (dm_info), intent(in) :: dminfo
-! type (exchange_list), pointer :: cellsToSend, cellsToRecv
if(monotonicOn) then
-! call mpas_tracer_advection_mono_tend( tend, s_old, s_new, diag, grid, dt, dminfo, cellsToSend, cellsToRecv)
- call mpas_tracer_advection_mono_tend(tracers, uh, w, h, dt, zDistance, zWeightK, zWeightKm1, grid, tend_h, tend)
+ call mpas_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend)
else
- call mpas_tracer_advection_std_tend(tracers, uh, w, zDistance, zWeightK, zWeightKm1, grid, tend)
+ call mpas_tracer_advection_std_tend(tracers, uh, w, verticalCellSize, grid, tend)
endif
end subroutine mpas_tracer_advection_tend!}}}
Modified: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_mono.F 2012-03-08 23:21:56 UTC (rev 1609)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_mono.F 2012-03-09 03:20:50 UTC (rev 1610)
@@ -15,7 +15,7 @@
contains
- subroutine mpas_tracer_advection_mono_tend(tracers, uh, w, h, dt, zDistance, zWeightK, zWeightKm1, grid, tend_h, tend)!{{{
+ subroutine mpas_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -27,11 +27,9 @@
real (kind=RKIND), dimension(:,:), intent(in) :: uh !< Input: Thichness weighted velocitiy
real (kind=RKIND), dimension(:,:), intent(in) :: w !< Input: Vertical velocitiy
real (kind=RKIND), dimension(:,:), intent(in) :: h !< Input: Thickness
+ real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !< Input: Distance between vertical interfaces of a cell
real (kind=RKIND), dimension(:,:), intent(in) :: tend_h !< Input: Tendency for thickness field
real (kind=RKIND), intent(in) :: dt !< Input: Timestep
- real (kind=RKIND), dimension(:), intent(in) :: zDistance !< Input: Distance between vertical interfaces
- real (kind=RKIND), dimension(:), intent(in) :: zWeightK !< Input: Weight for tracers to map to edge, weight from K cell
- real (kind=RKIND), dimension(:), intent(in) :: zWeightKm1 !< Input: Weight for tracers to map to edge, weight from K-1 cell
type (mesh_type), intent(in) :: grid
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
@@ -43,6 +41,7 @@
real (kind=RKIND) :: coef_3rd_order, flux_upwind, tracer_min_new, tracer_max_new, tracer_upwind_new, scale_factor
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 :: adv_coefs, adv_coefs_3rd
real (kind=RKIND), dimension(:,:), pointer :: tracer_cur, tracer_new, upwind_tendency, inv_h_new, tracer_max, tracer_min
@@ -95,7 +94,6 @@
end do
end do
-
! 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
@@ -104,7 +102,7 @@
tracer_cur(k,iCell) = tracers(iTracer,k,iCell)
upwind_tendency(k, iCell) = 0.0
- !dwj 02/21/12 tracer_new is supposed to be the "new" tracer state. This allows bounds checks.
+ !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
end if
@@ -117,12 +115,13 @@
! Compute the high order vertical flux. Also determine bounds on tracer_cur.
do iCell = 1, nCells
k = 1
-! high_order_vert_flux(k,iCell) = 0.
tracer_max(k,iCell) = max(tracer_cur(k,iCell),tracer_cur(k+1,iCell))
tracer_min(k,iCell) = min(tracer_cur(k,iCell),tracer_cur(k+1,iCell))
k = 2
- high_order_vert_flux(k,iCell) = w(k,iCell)*(zWeightK(k)*tracer_cur(k,iCell)+zWeightKm1(k)*tracer_cur(k-1,iCell))
+ verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+ verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+ high_order_vert_flux(k,iCell) = w(k,iCell)*(verticalWeightK*tracer_cur(k,iCell)+verticalWeightKm1*tracer_cur(k-1,iCell))
tracer_max(k,iCell) = max(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
tracer_min(k,iCell) = min(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
@@ -135,12 +134,12 @@
end do
k = maxLevelCell(iCell)
- high_order_vert_flux(k,iCell) = w(k,iCell)*(zWeightK(k)*tracer_cur(k,iCell)+zWeightKm1(k)*tracer_cur(k-1,iCell))
+ verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+ verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
+ high_order_vert_flux(k,iCell) = w(k,iCell)*(verticalWeightK*tracer_cur(k,iCell)+verticalWeightKm1*tracer_cur(k-1,iCell))
tracer_max(k,iCell) = max(tracer_cur(k,iCell),tracer_cur(k-1,iCell))
tracer_min(k,iCell) = min(tracer_cur(k,iCell),tracer_cur(k-1,iCell))
-! high_order_vert_flux(maxLevelCell(iCell)+1,iCell) = 0.
-
! pull tracer_min and tracer_max from the (horizontal) surrounding cells
do i = 1, nEdgesOnCell(iCell)
do k=1, maxLevelCell(cellsOnCell(i, iCell))
@@ -167,7 +166,7 @@
! Upwind fluxes are accumulated in upwind_tendency
do iCell = 1, nCells
do k = 2, maxLevelCell(iCell)
- !dwj: wrong for ocean?
+ ! dwj 02/03/12 Ocean and Atmosphere are different in vertical
! flux_upwind = max(0.,w(k,iCell))*tracer_cur(k-1,iCell) + min(0.,w(k,iCell))*tracer_cur(k,iCell)
flux_upwind = min(0.,w(k,iCell))*tracer_cur(k-1,iCell) + max(0.,w(k,iCell))*tracer_cur(k,iCell)
upwind_tendency(k-1,iCell) = upwind_tendency(k-1,iCell) + flux_upwind
@@ -268,9 +267,11 @@
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
- !dwj 02/21/12 tracer_new holds a tendency for now.
- tracer_new(k, cell1) = tracer_new(k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
- tracer_new(k, cell2) = tracer_new(k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
+ if (config_check_monotonicity) then
+ !tracer_new holds a tendency for now.
+ tracer_new(k, cell1) = tracer_new(k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
+ tracer_new(k, cell2) = tracer_new(k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
+ end if
end do ! k loop
end do ! iEdge loop
@@ -280,23 +281,23 @@
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)
if (config_check_monotonicity) then
- !dwj 02/21/12 tracer_new holds a tendency for now.
+ !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)
- !dwj 02/21/12 now, updated tracer_new to be the new state of the tracer.
+ !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)
end if
end do ! k loop
end do ! iCell loop
if (config_check_monotonicity) then
- !dwj 02/21/12 build min and max of cur and new tracer states to check bounds
+ !build min and max bounds on old and new tracer for check on monotonicity.
cur_min = minval(tracer_cur(:,1:nCellsSolve))
cur_max = maxval(tracer_cur(:,1:nCellsSolve))
new_min = minval(tracer_new(:,1:nCellsSolve))
new_max = maxval(tracer_new(:,1:nCellsSolve))
- !dwj 02/21/12 do bounds check
+ !check monotonicity
if(new_min < cur_min-eps) then
write(*,*) 'Minimum out of bounds on tracer ', iTracer, cur_min, new_min
end if
Modified: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std.F 2012-03-08 23:21:56 UTC (rev 1609)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std.F 2012-03-09 03:20:50 UTC (rev 1610)
@@ -18,7 +18,7 @@
contains
- subroutine mpas_tracer_advection_std_tend(tracers, uh, w, zDistance, zWeightK, zWeightKm1, grid, tend)!{{{
+ subroutine mpas_tracer_advection_std_tend(tracers, uh, w, verticalCellSize, grid, tend)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -30,14 +30,14 @@
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
real (kind=RKIND), dimension(:,:), intent(in) :: uh, w
- real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !< Input: Distance between vertical interfaces of a cell
type (mesh_type), intent(in) :: grid
call mpas_timer_start("tracer-hadv", .false.)
call mpas_tracer_advection_std_hadv_tend(tracers, uh, grid, tend)
call mpas_timer_stop("tracer-hadv")
call mpas_timer_start("tracer-vadv", .false.)
- call mpas_tracer_advection_std_vadv_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)
+ call mpas_tracer_advection_std_vadv_tend(tracers, w, verticalCellSize, grid, tend)
call mpas_timer_stop("tracer-vadv")
end subroutine mpas_tracer_advection_std_tend!}}}
Modified: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv.F 2012-03-08 23:21:56 UTC (rev 1609)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv.F 2012-03-09 03:20:50 UTC (rev 1610)
@@ -20,21 +20,20 @@
contains
- subroutine mpas_tracer_advection_std_vadv_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)!{{{
+ subroutine mpas_tracer_advection_std_vadv_tend(tracers, w, verticalCellSize, grid, tend)!{{{
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
- real (kind=RKIND), dimension(:,:), intent(in) :: w
- real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ real (kind=RKIND), dimension(:,:), intent(in) :: w, verticalCellSize
type (mesh_type), intent(in) :: grid
real (kind=RKIND) :: dt
if(order2) then
- call mpas_tracer_advection_std_vadv2_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)
+ call mpas_tracer_advection_std_vadv2_tend(tracers, w, verticalCellSize, grid, tend)
else if(order3) then
- call mpas_tracer_advection_std_vadv3_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)
+ call mpas_tracer_advection_std_vadv3_tend(tracers, w, verticalCellSize, grid, tend)
else if(order4) then
- call mpas_tracer_advection_std_vadv4_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)
+ call mpas_tracer_advection_std_vadv4_tend(tracers, w, verticalCellSize, grid, tend)
endif
end subroutine mpas_tracer_advection_std_vadv_tend!}}}
Modified: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv2.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv2.F 2012-03-08 23:21:56 UTC (rev 1609)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv2.F 2012-03-09 03:20:50 UTC (rev 1610)
@@ -15,7 +15,7 @@
contains
- subroutine mpas_tracer_advection_std_vadv2_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)!{{{
+ subroutine mpas_tracer_advection_std_vadv2_tend(tracers, w, verticalCellSize, grid, tend)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -26,8 +26,7 @@
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
- real (kind=RKIND), dimension(:,:), intent(in) :: w
- real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ real (kind=RKIND), dimension(:,:), intent(in) :: w, verticalCellSize
type (mesh_type), intent(in) :: grid
integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
@@ -35,8 +34,8 @@
real (kind=RKIND) :: tracer_weight_cell1, tracer_weight_cell2
-! real (kind=RKIND), dimension( s_old % num_tracers, grid % nVertLevels + 1 ) :: vert_flux
real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
+ real (kind=RKIND) :: weightK, weightKm1
integer :: nVertLevels, num_tracers
integer, dimension(:), pointer :: maxLevelCell
@@ -59,7 +58,9 @@
do iCell=1,grid % nCellsSolve
do k = 2, maxLevelCell(iCell)
do iTracer=1,num_tracers
- vert_flux(iTracer,k) = w(k,iCell)*(zWeightK(k)*tracers(iTracer,k,iCell)+zWeightKm1(k)*tracers(iTracer,k-1,iCell))
+ weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
+ weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
+ vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*tracers(iTracer,k-1,iCell))
end do
end do
@@ -68,8 +69,6 @@
do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + ( vert_flux(iTracer, k+1) - vert_flux(iTracer, k))
-! tracers(iTracer,k,iCell) = ( tracers(iTracer,k,iCell)*h_old(k,iCell) &
-! + dt*( tracer_tend(iTracer,k,iCell) -rdnw(k)*(vert_flux(iTracer,k+1)-vert_flux(iTracer,k)) ) )/h_new(k,iCell)
end do
end do
end do
Modified: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv3.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv3.F 2012-03-08 23:21:56 UTC (rev 1609)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv3.F 2012-03-09 03:20:50 UTC (rev 1610)
@@ -17,7 +17,7 @@
contains
- subroutine mpas_tracer_advection_std_vadv3_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)!{{{
+ subroutine mpas_tracer_advection_std_vadv3_tend(tracers, w, verticalCellSize, grid, tend)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -28,15 +28,14 @@
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
- real (kind=RKIND), dimension(:,:), intent(in) :: w
- real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ real (kind=RKIND), dimension(:,:), intent(in) :: w, verticalCellSize
type (mesh_type), intent(in) :: grid
real (kind=RKIND) :: dt
integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
-! real (kind=RKIND), dimension( s_old % num_tracers, grid % nVertLevels + 1 ) :: vert_flux
real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
+ real (kind=RKIND) :: weightK, weightKm1
integer :: nVertLevels, num_tracers
integer, dimension(:), pointer :: maxLevelCell
@@ -51,13 +50,14 @@
allocate(vert_flux(num_tracers, nVertLevels+1))
vert_flux(:,1) = 0.
-! vert_flux(:,grid % nVertLevels+1) = 0.
do iCell=1,grid % nCellsSolve
k = 2
do iTracer=1,num_tracers
- vert_flux(iTracer,k) = w(k,iCell)*(zWeightK(k)*tracers(iTracer,k,iCell)+zWeightKm1(k)*tracers(iTracer,k-1,iCell))
+ weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
+ weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
+ vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*tracers(iTracer,k-1,iCell))
enddo
do k=3,maxLevelCell(iCell)-1
@@ -71,7 +71,9 @@
k = maxLevelCell(iCell)
do iTracer=1,num_tracers
- vert_flux(iTracer,k) = w(k,iCell)*(zWeightK(k)*tracers(iTracer,k,iCell)+zWeightKm1(k)*tracers(iTracer,k-1,iCell))
+ weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
+ weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
+ vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*tracers(iTracer,k-1,iCell))
enddo
vert_Flux(:, maxLevelCell(iCell)+1) = 0.0
@@ -79,9 +81,6 @@
do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (vert_flux(iTracer, k+1) - vert_flux(iTracer, k))
-
-! tracers(iTracer,k,iCell) = ( tracers(iTracer,k,iCell)*h_old(k,iCell) &
-! + dt*( tracer_tend(iTracer,k,iCell) -rdnw(k)*(vert_flux(iTracer,k+1)-vert_flux(iTracer,k)) ) )/h_new(k,iCell)
end do
end do
end do
Modified: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv4.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv4.F 2012-03-08 23:21:56 UTC (rev 1609)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv4.F 2012-03-09 03:20:50 UTC (rev 1610)
@@ -15,7 +15,7 @@
contains
- subroutine mpas_tracer_advection_std_vadv4_tend(tracers, w, zDistance, zWeightK, zWeightKm1, grid, tend)!{{{
+ subroutine mpas_tracer_advection_std_vadv4_tend(tracers, w, verticalCellSize, grid, tend)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -26,15 +26,14 @@
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
- real (kind=RKIND), dimension(:,:), intent(in) :: w
- real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ real (kind=RKIND), dimension(:,:), intent(in) :: w, verticalCellSize
type (mesh_type), intent(in) :: grid
real (kind=RKIND) :: dt
integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
-! real (kind=RKIND), dimension( s_old % num_tracers, grid % nVertLevels + 1 ) :: vert_flux
real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
+ real (kind=RKIND) :: weightK, weightKm1
integer :: nVertLevels, num_tracers
integer, dimension(:), pointer :: maxLevelCell
@@ -55,7 +54,9 @@
k = 2
do iTracer=1,num_tracers
- vert_flux(iTracer,k) = w(k,iCell)*(zWeightK(k)*tracers(iTracer,k,iCell)+zWeightKm1(k)*tracers(iTracer,k-1,iCell))
+ weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
+ weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
+ vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*tracers(iTracer,k-1,iCell))
enddo
do k=3,nVertLevels-1
do iTracer=1,num_tracers
@@ -66,7 +67,9 @@
k = maxLevelCell(iCell)
do iTracer=1,num_tracers
- vert_flux(iTracer,k) = w(k,iCell)*(zWeightK(k)*tracers(iTracer,k,iCell)+zWeightKm1(k)*tracers(iTracer,k-1,iCell))
+ weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
+ weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
+ vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*tracers(iTracer,k-1,iCell))
enddo
vert_flux(:,maxLevelCell(iCell)+1) = 0.0
@@ -74,8 +77,6 @@
do k=1,maxLevelCell(iCell)
do iTracer=1,num_tracers
tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (vert_flux(iTracer, k+1) - vert_flux(iTracer, k))
-! tracers(iTracer,k,iCell) = ( tracers(iTracer,k,iCell)*h_old(k,iCell) &
-! + dt*( tracer_tend(iTracer,k,iCell) -rdnw(k)*(vert_flux(iTracer,k+1)-vert_flux(iTracer,k)) ) )/h_new(k,iCell)
end do
end do
</font>
</pre>