<p><b>dwj07@fsu.edu</b> 2012-02-03 11:23:05 -0700 (Fri, 03 Feb 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Working version of monotonic transport. At least in horizontal. Vertical needs to be tested.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/advection_routines/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/advection_routines/src/core_ocean/Makefile 2012-02-03 18:14:22 UTC (rev 1461)
+++ branches/ocean_projects/advection_routines/src/core_ocean/Makefile 2012-02-03 18:23:05 UTC (rev 1462)
@@ -136,8 +136,7 @@
mpas_ocn_equation_of_state_linear.o:
-mpas_ocn_mpas_core.o: mpas_ocn_mpas_core.o \
- mpas_ocn_test_cases.o \
+mpas_ocn_mpas_core.o: mpas_ocn_test_cases.o \
mpas_ocn_advection.o \
mpas_ocn_thick_hadv.o \
mpas_ocn_thick_vadv.o \
Modified: branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_tendency.F 2012-02-03 18:14:22 UTC (rev 1461)
+++ branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_tendency.F 2012-02-03 18:23:05 UTC (rev 1462)
@@ -291,7 +291,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tend_scalar(tend, s, d, grid)!{{{
+ subroutine ocn_tend_scalar(tend, s, d, grid, dt, parinfo, dminfo)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -309,10 +309,13 @@
! type (diagnostics_type), intent(in) :: d ! dwj: 01/17/12 change for cross core compatibilty
type (diag_type), intent(in) :: d
type (mesh_type), intent(in) :: grid
+ real (kind=RKIND), intent(in) :: dt
+ type (parallel_info), intent(in) :: parinfo
+ type (dm_info), intent(in) :: dminfo
real (kind=RKIND), dimension(:), pointer :: hRatioZLevelK, hRatioZLevelKm1, zTopZLevel, hMeanTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
- u,h,wTop, h_edge, vertDiffTopOfCell
+ u,h,wTop, h_edge, vertDiffTopOfCell, tend_h
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
@@ -335,6 +338,7 @@
hMeanTopZLevel => grid % hMeanTopZLevel % array
tend_tr => tend % tracers % array
+ tend_h => tend % h % array
allocate(uh(grid % nVertLevels, grid % nEdges+1))
@@ -357,21 +361,21 @@
! 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)
-! call mpas_tracer_advection_tend(tracers, uh, wTop, hMeanTopZLevel, hRatioZLevelKm1, hRatioZLevelK, grid, tend_tr)
-! call mpas_timer_stop("adv", tracerHadvTimer)
- call mpas_timer_start("hadv", .false., tracerHadvTimer)
- call ocn_tracer_hadv_tend(grid, u, h_edge, tracers, tend_tr, err)
- call mpas_timer_stop("hadv", tracerHadvTimer)
+ call mpas_timer_start("adv", .false., tracerHadvTimer)
+ call mpas_tracer_advection_tend(tracers, uh, wTop, h, dt, hMeanTopZLevel, hRatioZLevelKm1, hRatioZLevelK, grid, tend_h, parinfo, dminfo, tend_tr)
+ call mpas_timer_stop("adv", tracerHadvTimer)
+! call mpas_timer_start("hadv", .false., tracerHadvTimer)
+! call ocn_tracer_hadv_tend(grid, u, h_edge, tracers, tend_tr, err)
+! call mpas_timer_stop("hadv", tracerHadvTimer)
!
! tracer tendency: vertical advection term -d/dz( h \phi w)
!
- call mpas_timer_start("vadv", .false., tracerVadvTimer)
- call ocn_tracer_vadv_tend(grid, wtop, tracers, tend_tr, err)
- call mpas_timer_stop("vadv", tracerVadvTimer)
+! call mpas_timer_start("vadv", .false., tracerVadvTimer)
+! call ocn_tracer_vadv_tend(grid, wtop, tracers, tend_tr, err)
+! call mpas_timer_stop("vadv", tracerVadvTimer)
!
! tracer tendency: del2 horizontal tracer diffusion, div(h \kappa_2 </font>
<font color="gray">abla \phi)
@@ -536,9 +540,9 @@
h_edge = -1e34
-! 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
+ 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,nEdges*hadv2nd
Modified: branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-02-03 18:14:22 UTC (rev 1461)
+++ branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-02-03 18:23:05 UTC (rev 1462)
@@ -179,7 +179,8 @@
call filter_btr_mode_tend_u(block % tend, provis, block % diag, block % mesh)
endif
- call ocn_tend_scalar(block % tend, provis, block % diag, block % mesh)
+ !call ocn_tend_scalar(block % tend, provis, block % diag, block % mesh, rk_substep_weights(rk_step), block % parinfo, block % domain % dminfo, block % state % time_levs(1) % state % h % array)
+ call ocn_tend_scalar(block % tend, provis, block % diag, block % mesh, dt, block % parinfo, block % domain % dminfo)
block => block % next
end do
call mpas_timer_stop("RK4-tendency computations")
@@ -229,10 +230,6 @@
provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
- if (config_prescribe_thickness) then
- provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
- end if
-
call ocn_diagnostic_solve(dt, provis, block % mesh)
block => block % next
@@ -329,10 +326,6 @@
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
- if (config_prescribe_thickness) then
- block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
- end if
-
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
Modified: branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-03 18:14:22 UTC (rev 1461)
+++ branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-03 18:23:05 UTC (rev 1462)
@@ -777,7 +777,7 @@
call ocn_tend_h(block % tend, block % state % time_levs(2) % state , block % diag, block % mesh)
endif
- call ocn_tend_scalar(block % tend, block % state % time_levs(2) % state , block % diag, block % mesh)
+ call ocn_tend_scalar(block % tend, block % state % time_levs(2) % state , block % diag, block % mesh, dt, block % parinfo, block % domain % dminfo)
block => block % next
end do
@@ -1017,10 +1017,6 @@
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
- if(config_prescribe_thickness) then
- block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
- end if
-
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
Modified: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection.F 2012-02-03 18:14:22 UTC (rev 1461)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection.F 2012-02-03 18:23:05 UTC (rev 1462)
@@ -170,14 +170,18 @@
end subroutine mpas_tracer_advection_coefficients!}}}
-! subroutine mpas_tracer_advection_tend(tracers, uh, w, zDistance, zWeightUp, zWeightDown, grid, dt, dminfo, cellsToSend, cellsToRecv, tend)
- subroutine mpas_tracer_advection_tend(tracers, uh, w, zDistance, zWeightUp, zWeightDown, grid, tend)!{{{
+! 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, parinfo, dminfo, tend)!{{{
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, zWeightUp, zWeightDown
+ real (kind=RKIND), dimension(:,:), intent(in) :: uh, w, h
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
+ real (kind=RKIND), intent(in) :: dt
type (mesh_type), intent(in) :: grid
+ type (parallel_info), intent(in) :: parinfo
+ type (dm_info), intent(in) :: dminfo
+ real (kind=RKIND), dimension(:,:), intent(in) :: tend_h
! real (kind=RKIND) :: dt
! type (dm_info), intent(in) :: dminfo
@@ -185,9 +189,9 @@
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, zDistance, zWeightUp, zWeightDown, grid, tend)
+ call mpas_tracer_advection_mono_tend(tracers, uh, w, h, dt, zDistance, zWeightK, zWeightKm1, grid, tend_h, parinfo, dminfo, tend)
else
- call mpas_tracer_advection_std_tend(tracers, uh, w, zDistance, zWeightUp, zWeightDown, grid, tend)
+ call mpas_tracer_advection_std_tend(tracers, uh, w, zDistance, zWeightK, zWeightKm1, grid, tend)
endif
end subroutine mpas_tracer_advection_tend!}}}
Modified: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_mono.F 2012-02-03 18:14:22 UTC (rev 1461)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_mono.F 2012-02-03 18:23:05 UTC (rev 1462)
@@ -1,5 +1,6 @@
module mpas_tracer_advection_mono
+ use mpas_kind_types
use mpas_grid_types
use mpas_configure
use mpas_dmpar
@@ -15,7 +16,7 @@
contains
! subroutine mpas_tracer_advection_mono_tend( tend, s_old, s_new, diag, grid, dt, dminfo, cellsToSend, cellsToRecv)
- subroutine mpas_tracer_advection_mono_tend(tracers, uh, w, zDistance, zWeightUp, zWeightDown, grid, tend)!{{{
+ subroutine mpas_tracer_advection_mono_tend(tracers, uh, w, h, dt, zDistance, zWeightK, zWeightKm1, grid, tend_h, parinfo, dminfo, tend)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -24,15 +25,17 @@
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
- real (kind=RKIND), dimension(:,:), intent(in) :: uh, w
- real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightUp, zWeightDown
+ real (kind=RKIND), dimension(:,:), intent(in) :: uh, w, h, tend_h
+ real (kind=RKIND), intent(in) :: dt
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightK, zWeightKm1
type (mesh_type), intent(in) :: grid
+ type (parallel_info), intent(in) :: parinfo
+ type (dm_info), intent(in) :: dminfo
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
real (kind=RKIND) :: flux, tracer_weight
- real (kind=RKIND), dimension(:,:,:), pointer :: tracer_cur_in, tracer_next_in
! real (kind=RKIND), dimension(:,:), pointer :: h_old, h_new
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
integer, dimension(:,:), pointer :: cellsOnEdge
@@ -41,8 +44,8 @@
integer, dimension(:), pointer :: nAdvCellsForEdge
real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_3rd
- real (kind=RKiND), dimension( grid % nVertLevels, grid % nCells ) :: tracer_cur, tracer_next, tendency_temp
- real (kind=RKIND), dimension( grid % nVertLevels, grid % nCells ) :: s_max, s_min, s_update
+ real (kind=RKiND), dimension( grid % nVertLevels, grid % nCells ) :: tracer_cur, tracer_next, tendency_temp, h_new
+ real (kind=RKIND), dimension( grid % nVertLevels, grid % nCells ) :: s_max, s_min
real (kind=RKIND), dimension( grid % nVertLevels, grid % nCells ) :: scale_in, scale_out
real (kind=RKIND), dimension( grid % nVertLevels, grid % nEdges ) :: horiz_flux
@@ -61,7 +64,7 @@
integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop
integer, parameter :: hadv_opt = 2
- real (kind=RKIND), parameter :: eps=1.e-20
+ real (kind=RKIND), parameter :: eps = 1.e-20
logical, parameter :: debug_print = .false.
coef_3rd_order = config_coef_3rd_order
@@ -71,7 +74,7 @@
! h_old => s_old % rho_zz % array
areaCell => grid % areaCell % array
-! rdnw => grid % rdzw % array
+! rdnw => grid % rdzw % array
nAdvCellsForEdge => grid % nAdvCellsForEdge % array
advCellsForEdge => grid % advCellsForEdge % array
@@ -84,22 +87,23 @@
num_tracers = size(tracers,dim=1)
-
- !
- ! Runge Kutta integration, so we compute fluxes from tracer_next values, update starts from tracer_cur
- !
-
! do one tracer at a time
-
do iTracer = 1, num_tracers
do iCell = 1, grid%nCells
- do k=1, grid%nVertLevels
+ do k=1, maxLevelCell(iCell)
tracer_cur(k,iCell) = tracers(iTracer,k,iCell)
tracer_next(k,iCell) = tracers(iTracer,k,iCell)
tendency_temp(k, iCell) = 0.0
+ h_new(k, iCell) = h(k, iCell) + dt * tend_h(k, iCell)
+
+ s_min(k, iCell) = tracer_cur(k, iCell)
+ s_max(k, iCell) = tracer_cur(k, iCell)
end do
end do
+! vert_flux = 0.0
+! horiz_flux = 0.0
+
! scmin = tracer_cur(1,1)
! scmax = tracer_cur(1,1)
! do iCell = 1, grid%nCells
@@ -116,32 +120,31 @@
!
- do iCell=1,grid % nCellsSolve
+ do iCell=1,grid % nCells
! zero flux at top and bottom
- vert_flux(1,iCell) = 0.
-
k = 1
- s_max(k,iCell) = max(tracer_cur(1,iCell),tracer_cur(2,iCell))
- s_min(k,iCell) = min(tracer_cur(1,iCell),tracer_cur(2,iCell))
+ vert_flux(k,iCell) = 0.
+! s_max(k,iCell) = max(tracer_cur(k,iCell),tracer_cur(k+1,iCell))
+! s_min(k,iCell) = min(tracer_cur(k,iCell),tracer_cur(k+1,iCell))
- k = 2
- vert_flux(k,iCell) = w(k,iCell)*(zWeightUp(k)*tracer_cur(k,iCell)+zWeightDown(k)*tracer_cur(k-1,iCell))
- s_max(k,iCell) = max(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
- s_min(k,iCell) = min(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
-
- do k=3,maxLevelCell(iCell)-1
- vert_flux(k,iCell) = flux3( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell), &
- tracer_cur(k ,iCell),tracer_cur(k+1,iCell), &
- w(k,iCell), coef_3rd_order )
- s_max(k,iCell) = max(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
- s_min(k,iCell) = min(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
- end do
-
- k = maxLevelCell(iCell)
- vert_flux(k,iCell) = w(k,iCell)*(zWeightUp(k)*tracer_cur(k,iCell)+zWeightDown(k)*tracer_cur(k-1,iCell))
- s_max(k,iCell) = max(tracer_cur(k,iCell),tracer_cur(k-1,iCell))
- s_min(k,iCell) = min(tracer_cur(k,iCell),tracer_cur(k-1,iCell))
+! k = 2
+! vert_flux(k,iCell) = w(k,iCell)*(zWeightK(k)*tracer_cur(k,iCell)+zWeightKm1(k)*tracer_cur(k-1,iCell))
+! s_max(k,iCell) = max(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
+! s_min(k,iCell) = min(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
+!
+! do k=3,maxLevelCell(iCell)-1
+! vert_flux(k,iCell) = flux3( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell), &
+! tracer_cur(k ,iCell),tracer_cur(k+1,iCell), &
+! w(k,iCell), coef_3rd_order )
+! s_max(k,iCell) = max(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
+! s_min(k,iCell) = min(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
+! end do
+!
+! k = maxLevelCell(iCell)
+! vert_flux(k,iCell) = w(k,iCell)*(zWeightK(k)*tracer_cur(k,iCell)+zWeightKm1(k)*tracer_cur(k-1,iCell))
+! s_max(k,iCell) = max(tracer_cur(k,iCell),tracer_cur(k-1,iCell))
+! s_min(k,iCell) = min(tracer_cur(k,iCell),tracer_cur(k-1,iCell))
vert_flux(maxLevelCell(iCell)+1,iCell) = 0.
@@ -181,12 +184,14 @@
! vertical flux divergence for upwind update, we will put upwind update into tracer_next, and put factor of dt in fluxes
- do iCell = 1, grid % nCellsSolve
+ do iCell = 1, grid % nCells
do k = 2, maxLevelCell(iCell)
- flux_upwind = max(0.,w(k,iCell))*tracer_cur(k-1,iCell) + min(0.,w(k,iCell))*tracer_cur(k,iCell)
- tendency_temp(k-1,iCell) = tendency_temp(k-1,iCell) - flux_upwind
- tendency_temp(k ,iCell) = tendency_temp(k ,iCell) + flux_upwind
+ !dwj: wrong for ocean?
+! 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)
+ tendency_temp(k-1,iCell) = tendency_temp(k-1,iCell) + flux_upwind
+ tendency_temp(k ,iCell) = tendency_temp(k ,iCell) - flux_upwind
vert_flux(k,iCell) = vert_flux(k,iCell) - flux_upwind
end do
@@ -194,8 +199,11 @@
! contributions to the update: first the vertical flux component, then the horizontal
do k=1,maxLevelCell(iCell)
- scale_in (k,iCell) = -(min(0.,vert_flux(k+1,iCell))-max(0.,vert_flux(k,iCell)))
- scale_out(k,iCell) = -(max(0.,vert_flux(k+1,iCell))-min(0.,vert_flux(k,iCell)))
+! scale_in (k,iCell) = -(min(0.,vert_flux(k+1,iCell))-max(0.,vert_flux(k,iCell)))
+! scale_out(k,iCell) = -(max(0.,vert_flux(k+1,iCell))-min(0.,vert_flux(k,iCell)))
+
+ scale_in (k, iCell) = max(0.0, vert_flux(k+1, iCell)) - min(0.0, vert_flux(k, iCell))
+ scale_out(k, iCell) = min(0.0, vert_flux(k+1, iCell)) - max(0.0, vert_flux(k, iCell))
end do
end do
@@ -210,27 +218,29 @@
if (cell1 <= grid%nCells .or. cell2 <= grid%nCells) then ! only for owned cells
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))
- horiz_flux(k,iEdge) = horiz_flux(k,iEdge) - flux_upwind
tendency_temp(k,cell1) = tendency_temp(k,cell1) - flux_upwind / areaCell(cell1)
tendency_temp(k,cell2) = tendency_temp(k,cell2) + flux_upwind / areaCell(cell2)
+ horiz_flux(k,iEdge) = horiz_flux(k,iEdge) - flux_upwind
scale_out(k,cell1) = scale_out(k,cell1) - max(0.,horiz_flux(k,iEdge)) / areaCell(cell1)
scale_in (k,cell1) = scale_in (k,cell1) - min(0.,horiz_flux(k,iEdge)) / areaCell(cell1)
scale_out(k,cell2) = scale_out(k,cell2) + min(0.,horiz_flux(k,iEdge)) / areaCell(cell2)
scale_in (k,cell2) = scale_in (k,cell2) + max(0.,horiz_flux(k,iEdge)) / areaCell(cell2)
-
end do
end if
end do
! next, the limiter
-
- do iCell = 1, grid % nCellsSolve
+ do iCell = 1, grid % nCells
do k = 1, maxLevelCell(iCell)
- s_min_update = tendency_temp(k,iCell)+scale_out(k,iCell)
- s_max_update = tendency_temp(k,iCell)+scale_in (k,iCell)
- s_upwind = tendency_temp(k,iCell)
+! s_min_update = tendency_temp(k,iCell)+scale_out (k,iCell)
+! s_max_update = tendency_temp(k,iCell)+scale_in (k,iCell)
+! s_upwind = tendency_temp(k,iCell)
+ s_min_update = (tracer_cur(k,iCell)*h(k,iCell) + dt*(tendency_temp(k,iCell)+scale_out(k,iCell)))/h_new(k,iCell)
+ s_max_update = (tracer_cur(k,iCell)*h(k,iCell) + dt*(tendency_temp(k,iCell)+scale_in(k,iCell)))/h_new(k,iCell)
+ s_upwind = (tracer_cur(k,iCell)*h(k,iCell) + dt*tendency_temp(k,iCell))/h_new(k,iCell)
+ !dwj Need to fix scale_factor, because s_max and s_min don't mean what we need them to.
scale_factor = (s_max(k,iCell)-s_upwind)/(s_max_update-s_upwind+eps)
scale_in(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
@@ -239,27 +249,28 @@
end do
end do
+
!
! communicate scale factors here
!
-! call dmpar_exch_halo_field2dReal( dminfo, &
+! call mpas_dmpar_exch_halo_field2d_real( dminfo, &
! scale_in(:,:), &
! grid % nVertLevels, &
! grid % nCells, &
-! cellsToSend, cellsToRecv )
-! call dmpar_exch_halo_field2dReal( dminfo, &
+! parinfo % cellsToSend, parinfo % cellsToRecv )
+! call mpas_dmpar_exch_halo_field2d_real( dminfo, &
! scale_out(:,:), &
! grid % nVertLevels, &
! grid % nCells, &
-! cellsToSend, cellsToRecv )
+! parinfo % cellsToSend, parinfo % cellsToRecv )
!
! rescale the fluxes
!
do iEdge = 1, grid % nEdges
cell1 = grid % cellsOnEdge % array(1,iEdge)
cell2 = grid % cellsOnEdge % array(2,iEdge)
- if (cell1 <= grid%nCells .or. cell2 <= grid%nCells) then
- do k = 1, maxLevelCell(iCell)
+ if (cell1 <= grid%nCellsSolve .or. cell2 <= grid%nCellsSolve) then
+ do k = 1, maxLevelEdgeTop(iEdge)
flux = horiz_flux(k,iEdge)
flux = max(0.,flux) * min(scale_out(k,cell1), scale_in(k,cell2)) &
+ min(0.,flux) * min(scale_in(k,cell1), scale_out(k,cell2))
@@ -267,14 +278,16 @@
end do
end if
end do
-
+
! rescale the vertical flux
- do iCell=1,grid % nCells
+ do iCell=1,grid % nCellsSolve
do k = 2, maxLevelCell(iCell)
flux = vert_flux(k,iCell)
- flux = max(0.,flux) * min(scale_out(k-1,iCell), scale_in(k ,iCell)) &
- + min(0.,flux) * min(scale_out(k ,iCell), scale_in(k-1,iCell))
+! flux = max(0.,flux) * min(scale_out(k-1,iCell), scale_in(k ,iCell)) &
+! + min(0.,flux) * min(scale_out(k ,iCell), scale_in(k-1,iCell))
+ flux = max(0.,flux) * min(scale_out(k ,iCell), scale_in(k-1,iCell)) &
+ + min(0.,flux) * min(scale_out(k-1,iCell), scale_in(k ,iCell))
vert_flux(k,iCell) = flux
end do
end do
@@ -284,12 +297,12 @@
do iEdge=1,grid%nEdges
cell1 = cellsOnEdge(1,iEdge)
cell2 = cellsOnEdge(2,iEdge)
- if (cell1 <= grid%nCells .or. cell2 <= grid%nCells) then ! only for owned cells
+ if (cell1 <= grid%nCellsSolve .or. cell2 <= grid%nCellsSolve) then ! only for owned cells
do k=1,maxLevelEdgeTop(iEdge)
! tracer_next(k,cell1) = tracer_next(k,cell1) - horiz_flux(k,iEdge)/areaCell(cell1)
! tracer_next(k,cell2) = tracer_next(k,cell2) + horiz_flux(k,iEdge)/areaCell(cell2)
tend(iTracer, k, cell1) = tend(iTracer, k, cell1) - horiz_flux(k, iEdge)/areaCell(cell1)
- tend(iTracer, k, cell2) = tend(iTracer, k, cell2) - horiz_flux(k, iEdge)/areaCell(cell2)
+ tend(iTracer, k, cell2) = tend(iTracer, k, cell2) + horiz_flux(k, iEdge)/areaCell(cell2)
end do
end if
end do
@@ -298,31 +311,11 @@
do k=1,maxLevelCell(iCell)
! tracer_next(k,iCell) = ( tracer_next(k,iCell) &
! + (-rdnw(k)*(vert_flux(k+1,iCell)-vert_flux(k,iCell)) ) )/h_new(k,iCell)
- tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (vert_flux(k+1, iCell) - vert_flux(k, iCell))
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (vert_flux(k+1, iCell) - vert_flux(k, iCell)) + tendency_temp(k,iCell)
+! tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + tendency_temp(k,iCell)
end do
end do
-! if(debug_print) then
-
-! scmin = tracer_next(1,1)
-! scmax = tracer_next(1,1)
-! do iCell = 1, grid%nCellsSolve
-! do k=1, grid%nVertLevels
-! scmax = max(scmax,tracer_next(k,iCell))
-! scmin = min(scmin,tracer_next(k,iCell))
-! if(s_max(k,iCell) < tracer_next(k,iCell)) then
-! write(32,*) ' over - k,iCell,s_min,s_max,tracer_next ',k,iCell,s_min(k,iCell),s_max(k,iCell),tracer_next(k,iCell)
-! end if
-! if(s_min(k,iCell) > tracer_next(k,iCell)) then
-! write(32,*) ' under - k,iCell,s_min,s_max,tracer_next ',k,iCell,s_min(k,iCell),s_max(k,iCell),tracer_next(k,iCell)
-! end if
-! enddo
-! enddo
-! write(0,*) ' scmin, scmax new out ',scmin,scmax
-! write(0,*) ' icell_min, k_min ',icellmax, kmax
-
-! end if
-
! do iCell = 1, grid%nCells
! do k=1, grid%maxLevelCell(iCell)
! tracer_next_in(iTracer,k,iCell) = max(0.,tracer_next(k,iCell))
</font>
</pre>