<p><b>dwj07@fsu.edu</b> 2012-01-23 11:21:38 -0700 (Mon, 23 Jan 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
"Working" version of monotonic transport.<br>
<br>
Needs to be verified, but it compiles and runs.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/advection_routines/src/operators/Makefile
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/Makefile 2012-01-23 18:10:00 UTC (rev 1409)
+++ branches/ocean_projects/advection_routines/src/operators/Makefile 2012-01-23 18:21:38 UTC (rev 1410)
@@ -4,6 +4,7 @@
mpas_vector_reconstruction.o \
mpas_spline_interpolation.o \
mpas_tracer_advection.o \
+ mpas_tracer_advection_mono.o \
mpas_tracer_advection_std.o \
mpas_tracer_advection_std_hadv.o \
mpas_tracer_advection_std_vadv.o \
@@ -23,8 +24,10 @@
mpas_spline_interpolation:
-mpas_tracer_advection.o: mpas_tracer_advection_std.o
+mpas_tracer_advection.o: mpas_tracer_advection_std.o mpas_tracer_advection_mono.o
+mpas_tracer_advection_mono.o: mpas_tracer_advection_helpers.o
+
mpas_tracer_advection_std.o: mpas_tracer_advection_std_hadv.o mpas_tracer_advection_std_vadv.o
mpas_tracer_advection_std_hadv.o: mpas_tracer_advection_helpers.o
Modified: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection.F 2012-01-23 18:10:00 UTC (rev 1409)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection.F 2012-01-23 18:21:38 UTC (rev 1410)
@@ -5,7 +5,7 @@
use mpas_configure
use mpas_tracer_advection_std
-! use mpas_tracer_advection_mono
+ use mpas_tracer_advection_mono
implicit none
private
@@ -185,8 +185,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)
! 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, zWeightUp, zWeightDown, 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-01-23 18:10:00 UTC (rev 1409)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_mono.F 2012-01-23 18:21:38 UTC (rev 1410)
@@ -4,13 +4,18 @@
use mpas_configure
use mpas_dmpar
+ use mpas_tracer_advection_helpers
+
implicit none
private
save
+ public :: mpas_tracer_advection_mono_tend
+
contains
- subroutine mpas_tracer_advection_mono_tend( tend, s_old, s_new, diag, grid, dt, dminfo, cellsToSend, cellsToRecv)!{{{
+! 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)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -18,178 +23,136 @@
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- type (tend_type), intent(in) :: tend
- type (state_type), intent(in) :: s_old
- type (state_type), intent(inout) :: s_new
- type (diag_type), intent(in) :: diag
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: uh, w
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightUp, zWeightDown
type (mesh_type), intent(in) :: grid
- real (kind=RKIND) :: dt
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
- type (dm_info), intent(in) :: dminfo
- type (exchange_list), pointer :: cellsToSend, cellsToRecv
+ integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
+ real (kind=RKIND) :: flux, tracer_weight
-
- integer :: i, iCell, iEdge, k, iScalar, cell1, cell2
- real (kind=RKIND) :: flux, scalar_edge, d2fdx2_cell1, d2fdx2_cell2, scalar_weight
- real (kind=RKIND) :: scalar_weight_cell1, scalar_weight_cell2
-
- real (kind=RKIND), dimension(:,:,:), pointer :: scalar_old_in, scalar_new_in, scalar_tend
- real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- real (kind=RKIND), dimension(:,:), pointer :: uhAvg, h_old, h_new, wwAvg, rho_edge, rho_zz, zgrid, kdiff
- real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell, qv_init
+ 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
integer, dimension(:,:), pointer :: advCellsForEdge
integer, dimension(:), pointer :: nAdvCellsForEdge
real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_3rd
- real (kind=RKIND), dimension( grid % nVertLevels, grid % nCells ) :: scalar_old, scalar_new
+ 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 ) :: scale_in, scale_out
- real (kind=RKIND), dimension( grid % nVertLevels, grid % nEdges ) :: flux_arr
- real (kind=RKIND), dimension( grid % nVertLevels + 1, grid % nCells ) :: wdtn
+ real (kind=RKIND), dimension( grid % nVertLevels, grid % nEdges ) :: horiz_flux
+ real (kind=RKIND), dimension( grid % nVertLevels + 1, grid % nCells ) :: vert_flux
- integer :: nVertLevels, isc, num_scalars, icellmax, kmax
+ integer :: nVertLevels, num_tracers, icellmax, kmax
- real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw, meshScalingDel2, meshScalingDel4
+! real (kind=RKIND), dimension(:), pointer :: rdnw
real (kind=RKIND) :: coef_3rd_order
- real (kind=RKIND) :: h_theta_eddy_visc2, v_theta_eddy_visc2, scalar_turb_flux, z1,z2,z3,z4,zm,z0,zp
+ real (kind=RKIND) :: tracer_turb_flux, z1,z2,z3,z4,zm,z0,zp
- real (kind=RKIND) :: flux3, flux4, flux_upwind
- real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3, scmin,scmax
+ real (kind=RKIND) :: flux_upwind
real (kind=RKIND) :: s_min_update, s_max_update, s_upwind, scale_factor
+ integer, dimension(:), pointer :: maxLevelCell, maxLevelEdgeTop
+
integer, parameter :: hadv_opt = 2
real (kind=RKIND), parameter :: eps=1.e-20
logical, parameter :: debug_print = .false.
- flux4(q_im2, q_im1, q_i, q_ip1, ua) = &
- ua*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
-
- flux3(q_im2, q_im1, q_i, q_ip1, ua, coef3) = &
- flux4(q_im2, q_im1, q_i, q_ip1, ua) + &
- coef3*abs(ua)*((q_ip1 - q_im2)-3.*(q_i-q_im1))/12.0
-
coef_3rd_order = config_coef_3rd_order
- scalar_old_in => s_old % scalars % array
- scalar_new_in => s_new % scalars % array
- kdiff => diag % kdiff % array
- deriv_two => grid % deriv_two % array
- uhAvg => diag % ruAvg % array
dvEdge => grid % dvEdge % array
- dcEdge => grid % dcEdge % array
cellsOnEdge => grid % cellsOnEdge % array
- scalar_tend => tend % scalars % array
- h_old => s_old % rho_zz % array
- h_new => s_new % rho_zz % array
- wwAvg => diag % wwAvg % array
+! h_old => s_old % rho_zz % array
areaCell => grid % areaCell % array
- fnm => grid % fzm % array
- fnp => grid % fzp % array
- rdnw => grid % rdzw % array
- meshScalingDel2 => grid % meshScalingDel2 % array
- meshScalingDel4 => grid % meshScalingDel4 % array
+! rdnw => grid % rdzw % array
nAdvCellsForEdge => grid % nAdvCellsForEdge % array
advCellsForEdge => grid % advCellsForEdge % array
adv_coefs => grid % adv_coefs % array
adv_coefs_3rd => grid % adv_coefs_3rd % array
+ maxLevelCell => grid % maxLevelCell % array
+ maxLevelEdgeTop => grid % maxLevelEdgeTop % array
nVertLevels = grid % nVertLevels
- h_theta_eddy_visc2 = config_h_theta_eddy_visc2
- v_theta_eddy_visc2 = config_v_theta_eddy_visc2
- rho_edge => diag % rho_edge % array
- rho_zz => s_new % rho_zz % array
- qv_init => grid % qv_init % array
- zgrid => grid % zgrid % array
+ num_tracers = size(tracers,dim=1)
-#ifndef DO_PHYSICS
- scalar_tend = 0. ! testing purposes - we have no sources or sinks
-#endif
!
- ! Runge Kutta integration, so we compute fluxes from scalar_new values, update starts from scalar_old
+ ! Runge Kutta integration, so we compute fluxes from tracer_next values, update starts from tracer_cur
!
- ! do one scalar at a time
+ ! do one tracer at a time
- num_scalars = 1
+ do iTracer = 1, num_tracers
+ write(0,*) ' mono transport for tracer ',iTracer
- do iScalar = 1, s_old % num_scalars
- write(0,*) ' mono transport for scalar ',iScalar
-
do iCell = 1, grid%nCells
- do k=1, grid%nVertLevels
- scalar_old(k,iCell) = scalar_old_in(iScalar,k,iCell)
- scalar_new(k,iCell) = scalar_new_in(iScalar,k,iCell)
+ do k=1, grid%nVertLevels
+ tracer_cur(k,iCell) = tracers(iTracer,k,iCell)
+ tracer_next(k,iCell) = tracers(iTracer,k,iCell)
+ tendency_temp(k, iCell) = 0.0
+ end do
end do
- end do
- scmin = scalar_old(1,1)
- scmax = scalar_old(1,1)
- do iCell = 1, grid%nCells
- do k=1, grid%nVertLevels
- scmin = min(scmin,scalar_old(k,iCell))
- scmax = max(scmax,scalar_old(k,iCell))
- enddo
- enddo
- write(0,*) ' scmin, scmin old in ',scmin,scmax
+! 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
- scmin = scalar_new(1,1)
- scmax = scalar_new(1,1)
- do iCell = 1, grid%nCells
- do k=1, grid%nVertLevels
- scmin = min(scmin,scalar_new(k,iCell))
- scmax = max(scmax,scalar_new(k,iCell))
- enddo
- enddo
- write(0,*) ' scmin, scmin new in ',scmin,scmax
+ !
+ ! vertical flux divergence, and min and max bounds for flux limiter
+ !
- !
- ! vertical flux divergence, and min and max bounds for flux limiter
- !
-
do iCell=1,grid % nCellsSolve
! zero flux at top and bottom
- wdtn(1,iCell) = 0.
- wdtn(grid % nVertLevels+1,iCell) = 0.
+ vert_flux(1,iCell) = 0.
k = 1
- s_max(k,iCell) = max(scalar_old(1,iCell),scalar_old(2,iCell))
- s_min(k,iCell) = min(scalar_old(1,iCell),scalar_old(2,iCell))
+ 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))
k = 2
- wdtn(k,iCell) = wwAvg(k,iCell)*(fnm(k)*scalar_new(k,iCell)+fnp(k)*scalar_new(k-1,iCell))
- s_max(k,iCell) = max(scalar_old(k-1,iCell),scalar_old(k,iCell),scalar_old(k+1,iCell))
- s_min(k,iCell) = min(scalar_old(k-1,iCell),scalar_old(k,iCell),scalar_old(k+1,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-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,nVertLevels-1
- wdtn(k,iCell) = flux3( scalar_new(k-2,iCell),scalar_new(k-1,iCell), &
- scalar_new(k ,iCell),scalar_new(k+1,iCell), &
- wwAvg(k,iCell), coef_3rd_order )
- s_max(k,iCell) = max(scalar_old(k-1,iCell),scalar_old(k,iCell),scalar_old(k+1,iCell))
- s_min(k,iCell) = min(scalar_old(k-1,iCell),scalar_old(k,iCell),scalar_old(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 = nVertLevels
- wdtn(k,iCell) = wwAvg(k,iCell)*(fnm(k)*scalar_new(k,iCell)+fnp(k)*scalar_new(k-1,iCell))
- s_max(k,iCell) = max(scalar_old(k,iCell),scalar_old(k-1,iCell))
- s_min(k,iCell) = min(scalar_old(k,iCell),scalar_old(k-1,iCell))
+ 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))
+ 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, grid % nVertLevels
- s_max(k,iCell) = max(s_max(k,iCell),scalar_old(k, grid % CellsOnCell % array(i,iCell)))
- s_min(k,iCell) = min(s_min(k,iCell),scalar_old(k, grid % CellsOnCell % array(i,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
@@ -198,19 +161,19 @@
!
! horizontal flux divergence
- flux_arr(:,:) = 0.
+ horiz_flux(:,:) = 0.
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
+ 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,grid % nVertLevels
- scalar_weight = uhAvg(k,iEdge)*(adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uhAvg(k,iEdge))*adv_coefs_3rd(i,iEdge))
- flux_arr(k,iEdge) = flux_arr(k,iEdge) + scalar_weight* scalar_new(k,iCell)
+ 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
@@ -218,27 +181,23 @@
end do
-! vertical flux divergence for upwind update, we will put upwind update into scalar_new, and put factor of dt in fluxes
+! 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
- k = 1
- scalar_new(k,iCell) = scalar_old(k,iCell)*h_old(k,iCell)
-
- do k = 2, nVertLevels
- scalar_new(k,iCell) = scalar_old(k,iCell)*h_old(k,iCell)
- flux_upwind = dt*(max(0.,wwAvg(k,iCell))*scalar_old(k-1,iCell) + min(0.,wwAvg(k,iCell))*scalar_old(k,iCell))
- scalar_new(k-1,iCell) = scalar_new(k-1,iCell) - flux_upwind*rdnw(k-1)
- scalar_new(k ,iCell) = scalar_new(k ,iCell) + flux_upwind*rdnw(k)
- wdtn(k,iCell) = dt*wdtn(k,iCell) - flux_upwind
+ 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
+ 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
- do k=1,nVertLevels
- scale_in (k,iCell) = - rdnw(k)*(min(0.,wdtn(k+1,iCell))-max(0.,wdtn(k,iCell)))
- scale_out(k,iCell) = - rdnw(k)*(max(0.,wdtn(k+1,iCell))-min(0.,wdtn(k,iCell)))
+ 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)))
end do
end do
@@ -250,18 +209,17 @@
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,grid % nVertLevels
- flux_upwind = grid % dvEdge % array(iEdge) * dt * &
- (max(0.,uhAvg(k,iEdge))*scalar_old(k,cell1) + min(0.,uhAvg(k,iEdge))*scalar_old(k,cell2))
- flux_arr(k,iEdge) = dt*flux_arr(k,iEdge) - flux_upwind
- scalar_new(k,cell1) = scalar_new(k,cell1) - flux_upwind / areaCell(cell1)
- scalar_new(k,cell2) = scalar_new(k,cell2) + flux_upwind / areaCell(cell2)
+ 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)
- scale_out(k,cell1) = scale_out(k,cell1) - max(0.,flux_arr(k,iEdge)) / areaCell(cell1)
- scale_in (k,cell1) = scale_in (k,cell1) - min(0.,flux_arr(k,iEdge)) / areaCell(cell1)
- scale_out(k,cell2) = scale_out(k,cell2) + min(0.,flux_arr(k,iEdge)) / areaCell(cell2)
- scale_in (k,cell2) = scale_in (k,cell2) + max(0.,flux_arr(k,iEdge)) / areaCell(cell2)
+ 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
@@ -270,10 +228,10 @@
! next, the limiter
do iCell = 1, grid % nCellsSolve
- do k = 1, nVertLevels
- s_min_update = (scalar_new(k,iCell)+scale_out(k,iCell))/h_new(k,iCell)
- s_max_update = (scalar_new(k,iCell)+scale_in (k,iCell))/h_new(k,iCell)
- s_upwind = scalar_new(k,iCell)/h_new(k,iCell)
+ 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)
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) )
@@ -286,28 +244,28 @@
!
! communicate scale factors here
!
- call dmpar_exch_halo_field2dReal( dminfo, &
- scale_in(:,:), &
- grid % nVertLevels, &
- grid % nCells, &
- cellsToSend, cellsToRecv )
- call dmpar_exch_halo_field2dReal( dminfo, &
- scale_out(:,:), &
- grid % nVertLevels, &
- grid % nCells, &
- cellsToSend, cellsToRecv )
+! call dmpar_exch_halo_field2dReal( dminfo, &
+! scale_in(:,:), &
+! grid % nVertLevels, &
+! grid % nCells, &
+! cellsToSend, cellsToRecv )
+! call dmpar_exch_halo_field2dReal( dminfo, &
+! scale_out(:,:), &
+! grid % nVertLevels, &
+! grid % nCells, &
+! cellsToSend, 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, nVertLevels
- flux = flux_arr(k,iEdge)
+ if (cell1 <= grid%nCells .or. cell2 <= grid%nCells) then
+ do k = 1, maxLevelCell(iCell)
+ 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))
- flux_arr(k,iEdge) = flux
+ horiz_flux(k,iEdge) = flux
end do
end if
end do
@@ -315,62 +273,65 @@
! rescale the vertical flux
do iCell=1,grid % nCells
- do k = 2, nVertLevels
- flux = wdtn(k,iCell)
+ 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))
- wdtn(k,iCell) = flux
+ vert_flux(k,iCell) = flux
end do
end do
!
-! do the scalar update now that we have the fluxes
+! 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,grid % nVertLevels
- scalar_new(k,cell1) = scalar_new(k,cell1) - flux_arr(k,iEdge)/areaCell(cell1)
- scalar_new(k,cell2) = scalar_new(k,cell2) + flux_arr(k,iEdge)/areaCell(cell2)
+ if (cell1 <= grid%nCells .or. cell2 <= grid%nCells) 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
do iCell=1,grid % nCellsSolve
- do k=1,grid % nVertLevels
- scalar_new(k,iCell) = ( scalar_new(k,iCell) &
- + (-rdnw(k)*(wdtn(k+1,iCell)-wdtn(k,iCell)) ) )/h_new(k,iCell)
+ 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))
end do
end do
- if(debug_print) then
+! if(debug_print) then
- scmin = scalar_new(1,1)
- scmax = scalar_new(1,1)
- do iCell = 1, grid%nCellsSolve
- do k=1, grid%nVertLevels
- scmax = max(scmax,scalar_new(k,iCell))
- scmin = min(scmin,scalar_new(k,iCell))
- if(s_max(k,iCell) < scalar_new(k,iCell)) then
- write(32,*) ' over - k,iCell,s_min,s_max,scalar_new ',k,iCell,s_min(k,iCell),s_max(k,iCell),scalar_new(k,iCell)
- end if
- if(s_min(k,iCell) > scalar_new(k,iCell)) then
- write(32,*) ' under - k,iCell,s_min,s_max,scalar_new ',k,iCell,s_min(k,iCell),s_max(k,iCell),scalar_new(k,iCell)
- end if
- enddo
- enddo
- write(0,*) ' scmin, scmax new out ',scmin,scmax
- write(0,*) ' icell_min, k_min ',icellmax, kmax
+! 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
+! end if
- do iCell = 1, grid%nCells
- do k=1, grid%nVertLevels
- scalar_new_in(iScalar,k,iCell) = max(0.,scalar_new(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 ! loop over scalars
+ end do ! loop over tracers
end subroutine mpas_tracer_advection_mono_tend!}}}
</font>
</pre>