<p><b>dwj07@fsu.edu</b> 2012-02-03 12:28:49 -0700 (Fri, 03 Feb 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
General code clean up on monotonic advection.<br>
Improving Indentation.<br>
<br>
Removing halo updates for scale factors, as they are not required.<br>
</p><hr noshade><pre><font color="gray">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:30:43 UTC (rev 1463)
+++ branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_tendency.F 2012-02-03 19:28:49 UTC (rev 1464)
@@ -291,7 +291,7 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tend_scalar(tend, s, d, grid, dt, parinfo, dminfo)!{{{
+ subroutine ocn_tend_scalar(tend, s, d, grid, dt)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -310,8 +310,6 @@
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 :: &
@@ -362,7 +360,7 @@
! 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, h, dt, hMeanTopZLevel, hRatioZLevelKm1, hRatioZLevelK, grid, tend_h, parinfo, dminfo, tend_tr)
+ call mpas_tracer_advection_tend(tracers, uh, wTop, h, dt, hMeanTopZLevel, hRatioZLevelKm1, hRatioZLevelK, grid, tend_h, 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)
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:30:43 UTC (rev 1463)
+++ branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-02-03 19:28:49 UTC (rev 1464)
@@ -179,8 +179,7 @@
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, 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)
+ call ocn_tend_scalar(block % tend, provis, block % diag, block % mesh, dt)
block => block % next
end do
call mpas_timer_stop("RK4-tendency computations")
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:30:43 UTC (rev 1463)
+++ branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_time_integration_split.F 2012-02-03 19:28:49 UTC (rev 1464)
@@ -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, dt, block % parinfo, block % domain % dminfo)
+ call ocn_tend_scalar(block % tend, block % state % time_levs(2) % state , block % diag, block % mesh, dt)
block => block % next
end do
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:30:43 UTC (rev 1463)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection.F 2012-02-03 19:28:49 UTC (rev 1464)
@@ -171,7 +171,7 @@
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, parinfo, dminfo, tend)!{{{
+ subroutine mpas_tracer_advection_tend(tracers, uh, w, h, dt, zDistance, zWeightK, zWeightKm1, grid, tend_h, tend)!{{{
real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
@@ -179,8 +179,6 @@
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
@@ -189,7 +187,7 @@
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, parinfo, dminfo, tend)
+ call mpas_tracer_advection_mono_tend(tracers, uh, w, h, dt, zDistance, zWeightK, zWeightKm1, grid, tend_h, tend)
else
call mpas_tracer_advection_std_tend(tracers, uh, w, zDistance, zWeightK, zWeightKm1, grid, tend)
endif
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:30:43 UTC (rev 1463)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_mono.F 2012-02-03 19:28:49 UTC (rev 1464)
@@ -15,8 +15,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, h, dt, zDistance, zWeightK, zWeightKm1, grid, tend_h, parinfo, dminfo, tend)!{{{
+ subroutine mpas_tracer_advection_mono_tend(tracers, uh, w, h, dt, zDistance, zWeightK, zWeightKm1, grid, tend_h, tend)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -29,14 +28,11 @@
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 :: h_old, h_new
real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
integer, dimension(:,:), pointer :: cellsOnEdge
@@ -44,7 +40,7 @@
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, h_new
+ real (kind=RKiND), dimension( grid % nVertLevels, grid % nCells ) :: tracer_cur, 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
@@ -53,7 +49,6 @@
integer :: nVertLevels, num_tracers, icellmax, kmax
-! real (kind=RKIND), dimension(:), pointer :: rdnw
real (kind=RKIND) :: coef_3rd_order
real (kind=RKIND) :: tracer_turb_flux, z1,z2,z3,z4,zm,z0,zp
@@ -63,19 +58,14 @@
integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop
- integer, parameter :: hadv_opt = 2
real (kind=RKIND), parameter :: eps = 1.e-20
- logical, parameter :: debug_print = .false.
coef_3rd_order = config_coef_3rd_order
dvEdge => grid % dvEdge % array
cellsOnEdge => grid % cellsOnEdge % array
-! h_old => s_old % rho_zz % array
areaCell => grid % areaCell % array
-! rdnw => grid % rdzw % array
-
nAdvCellsForEdge => grid % nAdvCellsForEdge % array
advCellsForEdge => grid % advCellsForEdge % array
adv_coefs => grid % adv_coefs % array
@@ -87,12 +77,11 @@
num_tracers = size(tracers,dim=1)
- ! do one tracer at a time
+ ! 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
do iCell = 1, grid%nCells
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)
@@ -101,229 +90,181 @@
end do
end do
-! vert_flux = 0.0
-! horiz_flux = 0.0
+ vert_flux = 0.0
+ horiz_flux = 0.0
-! scmin = tracer_cur(1,1)
-! scmax = tracer_cur(1,1)
-! do iCell = 1, grid%nCells
-! do k=1, grid%nVertLevels
-! scmin = min(scmin,tracer_cur(k,iCell))
-! scmax = max(scmax,tracer_cur(k,iCell))
-! enddo
-! enddo
-! write(0,*) ' scmin, scmin old in ',scmin,scmax
-
-
!
- ! vertical flux divergence, and min and max bounds for flux limiter
+ ! Compute high order vertical flux. Also determine bounds on tracer_cur.
!
+ do iCell=1,grid % nCells
+ k = 1
+ 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))
-
- do iCell=1,grid % nCells
-
- ! zero flux at top and bottom
- k = 1
- 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)*(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))
+ 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
+ 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))
+ 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.
+ vert_flux(maxLevelCell(iCell)+1,iCell) = 0.
- ! pull s_min and s_max from the (horizontal) surrounding cells
-
- do i=1, grid % nEdgesOnCell % array(iCell)
- do k=1, maxLevelCell( grid % cellsOnCell % array(i, iCell))
- s_max(k,iCell) = max(s_max(k,iCell),tracer_cur(k, grid % CellsOnCell % array(i,iCell)))
- s_min(k,iCell) = min(s_min(k,iCell),tracer_cur(k, grid % CellsOnCell % array(i,iCell)))
- end do
- end do
-
- end do
-
- !
- ! horizontal flux divergence
-
- horiz_flux(:,:) = 0.
- 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
-
- do i=1,nAdvCellsForEdge(iEdge)
- iCell = advCellsForEdge(i,iEdge)
- do k=1,maxLevelCell(iCell)
- tracer_weight = uh(k,iEdge)*(adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
- horiz_flux(k,iEdge) = horiz_flux(k,iEdge) + tracer_weight* tracer_cur(k,iCell)
- end do
- end do
-
- end if
-
+ ! pull s_min and s_max from the (horizontal) surrounding cells
+ do i=1, grid % nEdgesOnCell % array(iCell)
+ do k=1, maxLevelCell( grid % cellsOnCell % array(i, iCell))
+ s_max(k,iCell) = max(s_max(k,iCell),tracer_cur(k, grid % CellsOnCell % array(i,iCell)))
+ s_min(k,iCell) = min(s_min(k,iCell),tracer_cur(k, grid % CellsOnCell % array(i,iCell)))
+ end do
end do
+ end do
-! 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 % nCells
-
- do k = 2, maxLevelCell(iCell)
- !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
-
-! scale_in(:,:) and scale_out(:,:) are used here to store the incoming and outgoing perturbation flux
-! contributions to the update: first the vertical flux component, then the horizontal
-
+ !
+ ! high order horizontal flux
+ !
+ horiz_flux(:,:) = 0.
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do i=1,nAdvCellsForEdge(iEdge)
+ iCell = advCellsForEdge(i,iEdge)
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) = 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))
+ tracer_weight = uh(k,iEdge)*(adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
+ horiz_flux(k,iEdge) = horiz_flux(k,iEdge) + tracer_weight* tracer_cur(k,iCell)
end do
+ end do
+ end do
+ !
+ ! low order upwind vertical flux (monotonic and diffused)
+ ! Remove low order flux from the high order flux.
+ ! Store left over high order flux in vert_flux array.
+ ! Upwind fluxes are accumulated in tendency_temp
+ !
+ do iCell = 1, grid % nCells
+ do k = 2, maxLevelCell(iCell)
+ !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
-! horizontal flux divergence for upwind update
+ ! scale_in contains the total remaining high order flux into iCell
+ ! it is positive.
+ ! scale_out contains the total remaining high order flux out of iCell
+ ! it is negative
+ do k=1,maxLevelCell(iCell)
+ ! dwj 02/03/12 Ocean and Atmosphere are different in vertical
+! 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)))
- ! upwind flux computation
+ 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
- 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
- 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))
- 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
+ !
+ ! low order upwind horizontal flux (monotinc and diffused)
+ ! Remove low order flux from the high order flux
+ ! Store left over high order flux in horiz_flux array
+ ! Upwind fluxes are accumulated in tendency_temp
+ !
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ 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))
+ 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
+ ! Accumulate remaining high order fluxes
+ 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 do
-! next, the limiter
- 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 = (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)
+ ! Build the factors for the FCT
+ do iCell = 1, grid % nCells
+ do k = 1, maxLevelCell(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) )
+ 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) )
- scale_factor = (s_upwind-s_min(k,iCell))/(s_upwind-s_min_update+eps)
- scale_out(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
+ scale_factor = (s_upwind-s_min(k,iCell))/(s_upwind-s_min_update+eps)
+ scale_out(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
- end do
end do
+ end do
-!
-! communicate scale factors here
-!
-! call mpas_dmpar_exch_halo_field2d_real( dminfo, &
-! scale_in(:,:), &
-! grid % nVertLevels, &
-! grid % nCells, &
-! parinfo % cellsToSend, parinfo % cellsToRecv )
-! call mpas_dmpar_exch_halo_field2d_real( dminfo, &
-! scale_out(:,:), &
-! grid % nVertLevels, &
-! grid % nCells, &
-! 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%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))
- horiz_flux(k,iEdge) = flux
- end do
- end if
- end do
+ ! rescale the high order horizontal fluxes
+ do iEdge = 1, grid % nEdges
+ cell1 = grid % cellsOnEdge % array(1,iEdge)
+ cell2 = grid % cellsOnEdge % array(2,iEdge)
+ 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))
+ horiz_flux(k,iEdge) = flux
+ end do
+ end do
- ! rescale the vertical flux
-
- 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 ,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
-!
-! do the tracer update now that we have the fluxes
-!
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- 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)
- end do
- end if
- end do
+ ! rescale the high order vertical flux
+ do iCell=1,grid % nCellsSolve
+ do k = 2, maxLevelCell(iCell)
+ flux = vert_flux(k,iCell)
+ ! dwj 02/03/12 Ocean and Atmosphere are different in vertical.
+! 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
- do iCell=1,grid % nCellsSolve
- 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)) + tendency_temp(k,iCell)
-! tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + tendency_temp(k,iCell)
- end do
+ ! Accumulate the scaled high order horizontal tendencies
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ do k=1,maxLevelEdgeTop(iEdge)
+ 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)
end do
+ end do
-! do iCell = 1, grid%nCells
+ ! Accumulate the scaled high order vertical tendencies, and the upwind tendencies
+ do iCell=1,grid % nCellsSolve
+ do k=1,maxLevelCell(iCell)
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (vert_flux(k+1, iCell) - vert_flux(k, iCell)) + tendency_temp(k,iCell)
+ end do
+ end do
+
+! do iCell = 1, grid%nCells
! do k=1, grid%maxLevelCell(iCell)
! tracer_next_in(iTracer,k,iCell) = max(0.,tracer_next(k,iCell))
! end do
-! end do
+! end do
end do ! loop over tracers
-
end subroutine mpas_tracer_advection_mono_tend!}}}
end module mpas_tracer_advection_mono
</font>
</pre>