<p><b>dwj07@fsu.edu</b> 2012-01-17 11:10:30 -0700 (Tue, 17 Jan 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Adding initial setup of tracer advection routines. They don't work yet.<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection.F (rev 0)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection.F 2012-01-17 18:10:30 UTC (rev 1382)
@@ -0,0 +1,212 @@
+module mpas_tracer_advection
+
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_constants
+ use mpas_dmpar
+ use mpas_vector_reconstruction
+ ! Added only clause to keep xlf90 from getting confused from the overloaded abs intrinsic in mpas_timekeeping
+ use mpas_mpas_timekeeping, only: MPAS_Time_type, MPAS_TimeInterval_type, &
+ MPAS_setTime, MPAS_setTimeInterval, MPAS_getTime, operator(+)
+
+ use mpas_tracer_advection_std
+ use mpas_tracer_advection_mono
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_init, &
+ mpas_tracer_advection_coefficients, &
+ mpas_tracer_advection_tend
+
+ logical :: monotonicOn
+
+ contains
+
+ subroutine mpas_tracer_advection_coefficients( grid )!{{{
+
+ implicit none
+ type (mesh_type) :: grid
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_3rd
+ integer, dimension(:,:), pointer :: cellsOnCell, cellsOnEdge, advCellsForEdge
+ integer, dimension(:), pointer :: nEdgesOnCell, nAdvCellsForEdge
+
+ integer :: cell1, cell2, iEdge, n, i, j, j_in, iCell
+ integer :: cell_list(20), ordered_cell_list(20)
+ logical :: addcell
+
+ deriv_two => grid % deriv_two % array
+ adv_coefs => grid % adv_coefs % array
+ adv_coefs_3rd => grid % adv_coefs_3rd % array
+ cellsOnCell => grid % cellsOnCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ advCellsForEdge => grid % advCellsForEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ nAdvCellsForEdge => grid % nAdvCellsForEdge % array
+
+ do iEdge = 1, grid % nEdges
+ nAdvCellsForEdge(iEdge) = 0
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ !
+ ! do only if this edge flux is needed to update owned cells
+ !
+ if (cell1 <= grid%nCells .or. cell2 <= grid%nCells) then
+
+ cell_list(1) = cell1
+ cell_list(2) = cell2
+ n = 2
+
+ ! add cells surrounding cell 1. n is number of cells currently in list
+ do i = 1, nEdgesOnCell(cell1)
+ if(cellsOnCell(i,cell1) /= cell2) then
+ n = n + 1
+ cell_list(n) = cellsOnCell(i,cell1)
+ end if
+ end do
+
+ ! add cells surrounding cell 2 (brute force approach)
+ do iCell = 1, nEdgesOnCell(cell2)
+ addcell = .true.
+ do i=1,n
+ if(cell_list(i) == cellsOnCell(iCell,cell2)) addcell = .false.
+ end do
+ if(addcell) then
+ n = n+1
+ cell_list(n) = cellsOnCell(iCell,cell2)
+ end if
+ end do
+
+ ! order the list by increasing cell number (brute force approach)
+
+ do i=1,n
+ ordered_cell_list(i) = grid % nCells + 2
+ j_in = 1
+ do j=1,n
+ if(ordered_cell_list(i) > cell_list(j) ) then
+ j_in = j
+ ordered_cell_list(i) = cell_list(j)
+ end if
+ end do
+! ordered_cell_list(i) = cell_list(j_in)
+ cell_list(j_in) = grid % nCells + 3
+ end do
+
+ nAdvCellsForEdge(iEdge) = n
+ do iCell = 1, nAdvCellsForEdge(iEdge)
+ advCellsForEdge(iCell,iEdge) = ordered_cell_list(iCell)
+ end do
+
+ ! we have the ordered list, now construct coefficients
+
+ adv_coefs(:,iEdge) = 0.
+ adv_coefs_3rd(:,iEdge) = 0.
+
+ ! pull together third and fourth order contributions to the flux
+ ! first from cell1
+
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cell1 ) j_in = j
+ end do
+ adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(1,1,iEdge)
+ adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) + deriv_two(1,1,iEdge)
+
+ do iCell = 1, nEdgesOnCell(cell1)
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cellsOnCell(iCell,cell1) ) j_in = j
+ end do
+ adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(iCell+1,1,iEdge)
+ adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) + deriv_two(iCell+1,1,iEdge)
+ end do
+
+ ! pull together third and fourth order contributions to the flux
+ ! now from cell2
+
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cell2 ) j_in = j
+ enddo
+ adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(1,2,iEdge)
+ adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) - deriv_two(1,2,iEdge)
+
+ do iCell = 1, nEdgesOnCell(cell2)
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cellsOnCell(iCell,cell2) ) j_in = j
+ enddo
+ adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(iCell+1,2,iEdge)
+ adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) - deriv_two(iCell+1,2,iEdge)
+ end do
+
+ do j = 1,n
+ adv_coefs (j,iEdge) = - (grid % dcEdge % array (iEdge) **2) * adv_coefs (j,iEdge) / 12.
+ adv_coefs_3rd(j,iEdge) = - (grid % dcEdge % array (iEdge) **2) * adv_coefs_3rd(j,iEdge) / 12.
+ end do
+
+ ! 2nd order centered contribution - place this in the main flux weights
+
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cell1 ) j_in = j
+ enddo
+ adv_coefs(j_in,iEdge) = adv_coefs(j_in,iEdge) + 0.5
+
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cell2 ) j_in = j
+ enddo
+ adv_coefs(j_in,iEdge) = adv_coefs(j_in,iEdge) + 0.5
+
+ ! multiply by edge length - thus the flux is just dt*ru times the results of the vector-vector multiply
+
+ do j=1,n
+ adv_coefs (j,iEdge) = grid % dvEdge % array(iEdge) * adv_coefs (j,iEdge)
+ adv_coefs_3rd(j,iEdge) = grid % dvEdge % array(iEdge) * adv_coefs_3rd(j,iEdge)
+ end do
+
+ end if ! only do for edges of owned-cells
+
+ end do ! end loop over edges
+
+ end subroutine adv_coef_compression!}}}
+
+ subroutine mpas_tracer_advection_tend(tend, s_old, s_new, diag, grid, dt, dminfo, cellsToSend, cellsToRecv)!{{{
+
+ 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
+ type (mesh_type), intent(in) :: grid
+ 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, s_old, s_new, diag, grid, dt, dminfo, cellsToSend, cellsToRecv)
+ else
+ call mpas_tracer_advection_std(tend, s_old, s_new, diag, grid, dt)
+ endif
+ end subroutine!}}}
+
+ subroutine mpas_tracer_advection_init(err)!{{{
+
+ integer, intent(inout) :: err
+
+ err = 0
+
+ monotonicOn = .false.
+
+ if(config_monotonic .and. config_positive_definite) then
+ monotonicOn = .true.
+ endif
+
+ end subroutine mpas_tracer_advection_init!}}}
+
+
+end module mpas_tracer_advection
Added: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_mono.F (rev 0)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_mono.F 2012-01-17 18:10:30 UTC (rev 1382)
@@ -0,0 +1,378 @@
+module mpas_tracer_advection_mono
+
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_constants
+ use mpas_dmpar
+ use mpas_vector_reconstruction
+ ! Added only clause to keep xlf90 from getting confused from the overloaded abs intrinsic in mpas_timekeeping
+ use mpas_mpas_timekeeping, only: MPAS_Time_type, MPAS_TimeInterval_type, &
+ MPAS_setTime, MPAS_setTimeInterval, MPAS_getTime, operator(+)
+
+ contains
+
+ subroutine mpas_tracer_advection_mono_tend( tend, s_old, s_new, diag, grid, dt, dminfo, cellsToSend, cellsToRecv)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ implicit none
+
+ 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
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND) :: dt
+
+ type (dm_info), intent(in) :: dminfo
+ type (exchange_list), pointer :: cellsToSend, cellsToRecv
+
+
+ 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
+ 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 ) :: 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
+
+ integer :: nVertLevels, isc, num_scalars, icellmax, kmax
+
+ real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw, meshScalingDel2, meshScalingDel4
+ 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) :: flux3, flux4, flux_upwind
+ real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3, scmin,scmax
+ real (kind=RKIND) :: s_min_update, s_max_update, s_upwind, scale_factor
+
+ 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
+ areaCell => grid % areaCell % array
+
+ fnm => grid % fzm % array
+ fnp => grid % fzp % array
+ rdnw => grid % rdzw % array
+ meshScalingDel2 => grid % meshScalingDel2 % array
+ meshScalingDel4 => grid % meshScalingDel4 % array
+
+ nAdvCellsForEdge => grid % nAdvCellsForEdge % array
+ advCellsForEdge => grid % advCellsForEdge % array
+ adv_coefs => grid % adv_coefs % array
+ adv_coefs_3rd => grid % adv_coefs_3rd % 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
+
+#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
+ !
+
+ ! do one scalar at a time
+
+ num_scalars = 1
+
+ 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)
+ 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 = 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
+ !
+
+
+ do iCell=1,grid % nCellsSolve
+
+ ! zero flux at top and bottom
+ wdtn(1,iCell) = 0.
+ wdtn(grid % nVertLevels+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))
+
+ 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))
+
+ 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))
+ 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))
+
+ ! 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)))
+ end do
+ end do
+
+ end do
+
+ !
+ ! horizontal flux divergence
+
+ flux_arr(:,:) = 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
+
+ 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)
+ end do
+ end do
+
+ end if
+
+ end do
+
+! vertical flux divergence for upwind update, we will put upwind update into scalar_new, 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
+ 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)))
+ end do
+
+ end do
+
+! horizontal flux divergence for upwind update
+
+ ! upwind flux computation
+
+ 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)
+
+ 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)
+
+ end do
+ end if
+ end do
+
+! 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)
+
+ 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) )
+
+ end do
+ end do
+!
+! 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 )
+!
+! 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)
+ 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
+ end do
+ end if
+ end do
+
+ ! rescale the vertical flux
+
+ do iCell=1,grid % nCells
+ do k = 2, nVertLevels
+ flux = wdtn(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
+ end do
+ end do
+!
+! do the scalar 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)
+ 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)
+ end do
+ end do
+
+ 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
+
+ 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
+
+ end do ! loop over scalars
+
+ end subroutine mpas_tracer_advection_mono_tend!}}}
+
+end module mpas_tracer_advection_mono
Added: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std.F (rev 0)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std.F 2012-01-17 18:10:30 UTC (rev 1382)
@@ -0,0 +1,385 @@
+module mpas_tracer_advection_std
+
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_constants
+ use mpas_dmpar
+ use mpas_vector_reconstruction
+ ! Added only clause to keep xlf90 from getting confused from the overloaded abs intrinsic in mpas_timekeeping
+ use mpas_mpas_timekeeping, only: MPAS_Time_type, MPAS_TimeInterval_type, &
+ MPAS_setTime, MPAS_setTimeInterval, MPAS_getTime, operator(+)
+
+ contains
+
+ subroutine mpas_tracer_advection_std_tend( tend, s_old, s_new, diag, grid, dt)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed scalar tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ implicit none
+
+ 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
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND) :: dt
+
+ 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, scalar_new, 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
+ 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( s_old % num_scalars, grid % nVertLevels ) :: flux_arr
+
+ real (kind=RKIND), dimension( s_old % num_scalars, grid % nVertLevels ) :: d2fdx2_cell1_arr, d2fdx2_cell2_arr
+
+ real (kind=RKIND), dimension( s_old % num_scalars, grid % nVertLevels + 1 ) :: wdtn
+ integer :: nVertLevels
+
+ real (kind=RKIND), dimension(:), pointer :: fnm, fnp, rdnw, meshScalingDel2, meshScalingDel4
+ 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) :: flux3, flux4
+ real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, ua, coef3
+
+ integer, parameter :: hadv_opt = 2
+
+ 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 => s_old % scalars % array
+ scalar_new => 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
+ areaCell => grid % areaCell % array
+
+ fnm => grid % fzm % array
+ fnp => grid % fzp % array
+ rdnw => grid % rdzw % array
+ meshScalingDel2 => grid % meshScalingDel2 % array
+ meshScalingDel4 => grid % meshScalingDel4 % array
+
+ nAdvCellsForEdge => grid % nAdvCellsForEdge % array
+ advCellsForEdge => grid % advCellsForEdge % array
+ adv_coefs => grid % adv_coefs % array
+ adv_coefs_3rd => grid % adv_coefs_3rd % 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
+
+#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 form scalar_old
+ !
+ !
+ ! horizontal flux divergence, accumulate in scalar_tend
+
+ 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
+
+ flux_arr(:,:) = 0.
+ do i=1,nAdvCellsForEdge(iEdge)
+ iCell = advCellsForEdge(i,iEdge)
+ do k=1,grid % nVertLevels
+ scalar_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uhAvg(k,iEdge))*adv_coefs_3rd(i,iEdge)
+ do iScalar=1,s_old % num_scalars
+ flux_arr(iScalar,k) = flux_arr(iScalar,k) + scalar_weight* scalar_new(iScalar,k,iCell)
+ end do
+ end do
+ end do
+
+ do k=1,grid % nVertLevels
+ do iScalar=1,s_old % num_scalars
+ scalar_tend(iScalar,k,cell1) = scalar_tend(iScalar,k,cell1) &
+ - uhAvg(k,iEdge)*flux_arr(iScalar,k)/areaCell(cell1)
+ scalar_tend(iScalar,k,cell2) = scalar_tend(iScalar,k,cell2) &
+ + uhAvg(k,iEdge)*flux_arr(iScalar,k)/areaCell(cell2)
+ end do
+ end do
+
+ end if
+ end do
+
+ !
+ ! vertical flux divergence
+ !
+
+ ! zero fluxes at top and bottom
+
+ wdtn(:,1) = 0.
+ wdtn(:,grid % nVertLevels+1) = 0.
+
+ if (config_scalar_vadv_order == 2) then
+
+ do iCell=1,grid % nCellsSolve
+ do k = 2, nVertLevels
+ do iScalar=1,s_old % num_scalars
+ wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
+ end do
+ end do
+ do k=1,grid % nVertLevelsSolve
+ do iScalar=1,s_old % num_scalars
+ scalar_new(iScalar,k,iCell) = ( scalar_old(iScalar,k,iCell)*h_old(k,iCell) &
+ + dt*( scalar_tend(iScalar,k,iCell) -rdnw(k)*(wdtn(iScalar,k+1)-wdtn(iScalar,k)) ) )/h_new(k,iCell)
+ end do
+ end do
+ end do
+
+ else if ( config_scalar_vadv_order == 3 ) then
+
+ do iCell=1,grid % nCellsSolve
+
+ k = 2
+ do iScalar=1,s_old % num_scalars
+ wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
+ enddo
+
+ do k=3,nVertLevels-1
+ do iScalar=1,s_old % num_scalars
+ wdtn(iScalar,k) = flux3( scalar_new(iScalar,k-2,iCell),scalar_new(iScalar,k-1,iCell), &
+ scalar_new(iScalar,k ,iCell),scalar_new(iScalar,k+1,iCell), &
+ wwAvg(k,iCell), coef_3rd_order )
+ end do
+ end do
+ k = nVertLevels
+ do iScalar=1,s_old % num_scalars
+ wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
+ enddo
+
+ do k=1,grid % nVertLevelsSolve
+ do iScalar=1,s_old % num_scalars
+ scalar_new(iScalar,k,iCell) = ( scalar_old(iScalar,k,iCell)*h_old(k,iCell) &
+ + dt*( scalar_tend(iScalar,k,iCell) -rdnw(k)*(wdtn(iScalar,k+1)-wdtn(iScalar,k)) ) )/h_new(k,iCell)
+ end do
+ end do
+
+ end do
+
+ else if ( config_scalar_vadv_order == 4 ) then
+
+ do iCell=1,grid % nCellsSolve
+
+ k = 2
+ do iScalar=1,s_old % num_scalars
+ wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
+ enddo
+ do k=3,nVertLevels-1
+ do iScalar=1,s_old % num_scalars
+ wdtn(iScalar,k) = flux4( scalar_new(iScalar,k-2,iCell),scalar_new(iScalar,k-1,iCell), &
+ scalar_new(iScalar,k ,iCell),scalar_new(iScalar,k+1,iCell), wwAvg(k,iCell) )
+ end do
+ end do
+ k = nVertLevels
+ do iScalar=1,s_old % num_scalars
+ wdtn(iScalar,k) = wwAvg(k,iCell)*(fnm(k)*scalar_new(iScalar,k,iCell)+fnp(k)*scalar_new(iScalar,k-1,iCell))
+ enddo
+
+ do k=1,grid % nVertLevelsSolve
+ do iScalar=1,s_old % num_scalars
+ scalar_new(iScalar,k,iCell) = ( scalar_old(iScalar,k,iCell)*h_old(k,iCell) &
+ + dt*( scalar_tend(iScalar,k,iCell) -rdnw(k)*(wdtn(iScalar,k+1)-wdtn(iScalar,k)) ) )/h_new(k,iCell)
+ end do
+ end do
+
+ end do
+
+ else
+
+ write(0,*) ' bad value for config_scalar_vadv_order - ', config_scalar_vadv_order
+
+ end if
+
+ end subroutine mpas_tracer_advection_std_tend!}}}
+
+ subroutine mpas_tracer_advection_coefficients( grid )!{{{
+
+ implicit none
+ type (mesh_type) :: grid
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_3rd
+ integer, dimension(:,:), pointer :: cellsOnCell, cellsOnEdge, advCellsForEdge
+ integer, dimension(:), pointer :: nEdgesOnCell, nAdvCellsForEdge
+
+ integer :: cell1, cell2, iEdge, n, i, j, j_in, iCell
+ integer :: cell_list(20), ordered_cell_list(20)
+ logical :: addcell
+
+ deriv_two => grid % deriv_two % array
+ adv_coefs => grid % adv_coefs % array
+ adv_coefs_3rd => grid % adv_coefs_3rd % array
+ cellsOnCell => grid % cellsOnCell % array
+ cellsOnEdge => grid % cellsOnEdge % array
+ advCellsForEdge => grid % advCellsForEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
+ nAdvCellsForEdge => grid % nAdvCellsForEdge % array
+
+ do iEdge = 1, grid % nEdges
+ nAdvCellsForEdge(iEdge) = 0
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ !
+ ! do only if this edge flux is needed to update owned cells
+ !
+ if (cell1 <= grid%nCells .or. cell2 <= grid%nCells) then
+
+ cell_list(1) = cell1
+ cell_list(2) = cell2
+ n = 2
+
+ ! add cells surrounding cell 1. n is number of cells currently in list
+ do i = 1, nEdgesOnCell(cell1)
+ if(cellsOnCell(i,cell1) /= cell2) then
+ n = n + 1
+ cell_list(n) = cellsOnCell(i,cell1)
+ end if
+ end do
+
+ ! add cells surrounding cell 2 (brute force approach)
+ do iCell = 1, nEdgesOnCell(cell2)
+ addcell = .true.
+ do i=1,n
+ if(cell_list(i) == cellsOnCell(iCell,cell2)) addcell = .false.
+ end do
+ if(addcell) then
+ n = n+1
+ cell_list(n) = cellsOnCell(iCell,cell2)
+ end if
+ end do
+
+ ! order the list by increasing cell number (brute force approach)
+
+ do i=1,n
+ ordered_cell_list(i) = grid % nCells + 2
+ j_in = 1
+ do j=1,n
+ if(ordered_cell_list(i) > cell_list(j) ) then
+ j_in = j
+ ordered_cell_list(i) = cell_list(j)
+ end if
+ end do
+! ordered_cell_list(i) = cell_list(j_in)
+ cell_list(j_in) = grid % nCells + 3
+ end do
+
+ nAdvCellsForEdge(iEdge) = n
+ do iCell = 1, nAdvCellsForEdge(iEdge)
+ advCellsForEdge(iCell,iEdge) = ordered_cell_list(iCell)
+ end do
+
+ ! we have the ordered list, now construct coefficients
+
+ adv_coefs(:,iEdge) = 0.
+ adv_coefs_3rd(:,iEdge) = 0.
+
+ ! pull together third and fourth order contributions to the flux
+ ! first from cell1
+
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cell1 ) j_in = j
+ end do
+ adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(1,1,iEdge)
+ adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) + deriv_two(1,1,iEdge)
+
+ do iCell = 1, nEdgesOnCell(cell1)
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cellsOnCell(iCell,cell1) ) j_in = j
+ end do
+ adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(iCell+1,1,iEdge)
+ adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) + deriv_two(iCell+1,1,iEdge)
+ end do
+
+ ! pull together third and fourth order contributions to the flux
+ ! now from cell2
+
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cell2 ) j_in = j
+ enddo
+ adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(1,2,iEdge)
+ adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) - deriv_two(1,2,iEdge)
+
+ do iCell = 1, nEdgesOnCell(cell2)
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cellsOnCell(iCell,cell2) ) j_in = j
+ enddo
+ adv_coefs (j_in,iEdge) = adv_coefs (j_in,iEdge) + deriv_two(iCell+1,2,iEdge)
+ adv_coefs_3rd(j_in,iEdge) = adv_coefs_3rd(j_in,iEdge) - deriv_two(iCell+1,2,iEdge)
+ end do
+
+ do j = 1,n
+ adv_coefs (j,iEdge) = - (grid % dcEdge % array (iEdge) **2) * adv_coefs (j,iEdge) / 12.
+ adv_coefs_3rd(j,iEdge) = - (grid % dcEdge % array (iEdge) **2) * adv_coefs_3rd(j,iEdge) / 12.
+ end do
+
+ ! 2nd order centered contribution - place this in the main flux weights
+
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cell1 ) j_in = j
+ enddo
+ adv_coefs(j_in,iEdge) = adv_coefs(j_in,iEdge) + 0.5
+
+ j_in = 0
+ do j=1, n
+ if( ordered_cell_list(j) == cell2 ) j_in = j
+ enddo
+ adv_coefs(j_in,iEdge) = adv_coefs(j_in,iEdge) + 0.5
+
+ ! multiply by edge length - thus the flux is just dt*ru times the results of the vector-vector multiply
+
+ do j=1,n
+ adv_coefs (j,iEdge) = grid % dvEdge % array(iEdge) * adv_coefs (j,iEdge)
+ adv_coefs_3rd(j,iEdge) = grid % dvEdge % array(iEdge) * adv_coefs_3rd(j,iEdge)
+ end do
+
+ end if ! only do for edges of owned-cells
+
+ end do ! end loop over edges
+
+ end subroutine adv_coef_compression!}}}
+
+end module mpas_tracer_advection_std
</font>
</pre>