<p><b>dwj07@fsu.edu</b> 2012-01-18 15:43:25 -0700 (Wed, 18 Jan 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Adding ported versions of the standard advection scheme.<br>
<br>
They are cleaned up and currently compile with the ocean core.<br>
<br>
They need to be checked for validity, and also at some point they need to be checked for cross compatibility with cores.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/advection_routines/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/advection_routines/src/core_ocean/Registry 2012-01-18 20:39:35 UTC (rev 1388)
+++ branches/ocean_projects/advection_routines/src/core_ocean/Registry 2012-01-18 22:43:25 UTC (rev 1389)
@@ -68,8 +68,10 @@
namelist character eos config_eos_type linear
namelist character advection config_vert_tracer_adv stencil
namelist integer advection config_vert_tracer_adv_order 4
+namelist integer advection config_tracer_vadv_order 2
namelist integer advection config_tracer_adv_order 2
namelist integer advection config_thickness_adv_order 2
+namelist real advection config_coef_3rd_order 1.0
namelist logical advection config_positive_definite false
namelist logical advection config_monotonic false
namelist logical restore config_restoreTS false
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-01-18 20:39:35 UTC (rev 1388)
+++ branches/ocean_projects/advection_routines/src/core_ocean/mpas_ocn_tendency.F 2012-01-18 22:43:25 UTC (rev 1389)
@@ -21,6 +21,8 @@
use mpas_constants
use mpas_timer
+ use mpas_tracer_advection
+
use ocn_thick_hadv
use ocn_thick_vadv
@@ -278,7 +280,7 @@
!***********************************************************************
!
-! routine ocn_tendSalar
+! routine ocn_tend_scalar
!
!> \brief Computes scalar tendency
!> \author Doug Jacobsen
@@ -296,7 +298,7 @@
! grid - grid metadata
! note: the variable s % tracers really contains the tracers,
! not tracers*h
- !
+ !
! Output: tend - computed scalar tendencies
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -308,11 +310,14 @@
type (diag_type), intent(in) :: d
type (mesh_type), intent(in) :: grid
+ real (kind=RKIND), dimension(:), pointer :: hRatioZLevelK, hRatioZLevelKm1, zTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
u,h,wTop, h_edge, vertDiffTopOfCell
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
+ real (kind=RKIND), dimension(:,:), allocatable :: uh
+
integer :: err
call mpas_timer_start("ocn_tend_scalar")
@@ -324,7 +329,15 @@
h_edge => s % h_edge % array
vertDiffTopOfCell => d % vertDiffTopOfCell % array
+ zTopZLevel => grid % zTopZLevel % array
+ hRatioZLevelK => grid % hRatioZLevelK % array
+ hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
+
tend_tr => tend % tracers % array
+
+ allocate(uh(grid % nEdges, grid % nVertLevels))
+
+ uh = u * h_edge
!
! initialize tracer tendency (RHS of tracer equation) to zero.
@@ -339,18 +352,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("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, zTopZlevel, hRatioZLevelk, hRatioZLevelKm1, 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)
!
! 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)
Modified: branches/ocean_projects/advection_routines/src/operators/Makefile
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/Makefile 2012-01-18 20:39:35 UTC (rev 1388)
+++ branches/ocean_projects/advection_routines/src/operators/Makefile 2012-01-18 22:43:25 UTC (rev 1389)
@@ -1,6 +1,16 @@
.SUFFIXES: .F .o
-OBJS = mpas_rbf_interpolation.o mpas_vector_reconstruction.o mpas_spline_interpolation.o mpas_tracer_advection.o
+OBJS = mpas_rbf_interpolation.o \
+ mpas_vector_reconstruction.o \
+ mpas_spline_interpolation.o \
+ mpas_tracer_advection.o \
+ mpas_tracer_advection_std.o \
+ mpas_tracer_advection_std_hadv.o \
+ mpas_tracer_advection_std_vadv.o \
+ mpas_tracer_advection_std_vadv2.o \
+ mpas_tracer_advection_std_vadv3.o \
+ mpas_tracer_advection_std_vadv4.o \
+ mpas_tracer_advection_helpers.o
all: operators
@@ -13,8 +23,22 @@
mpas_spline_interpolation:
-mpas_tracer_advection.o:
+mpas_tracer_advection.o: mpas_tracer_advection_std.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
+
+mpas_tracer_advection_std_vadv.o: mpas_tracer_advection_std_vadv2.o mpas_tracer_advection_std_vadv3.o mpas_tracer_advection_std_vadv4.o
+
+mpas_tracer_advection_std_vadv2.o: mpas_tracer_advection_helpers.o
+
+mpas_tracer_advection_std_vadv3.o: mpas_tracer_advection_helpers.o
+
+mpas_tracer_advection_std_vadv4.o: mpas_tracer_advection_helpers.o
+
+mpas_tracer_advection_helpers.o:
+
clean:
$(RM) *.o *.mod *.f90 libops.a
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-18 20:39:35 UTC (rev 1388)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection.F 2012-01-18 22:43:25 UTC (rev 1389)
@@ -1,9 +1,10 @@
module mpas_tracer_advection
+ use mpas_kind_types
use mpas_grid_types
use mpas_configure
-! use mpas_tracer_advection_std
+ use mpas_tracer_advection_std
! use mpas_tracer_advection_mono
implicit none
@@ -169,23 +170,24 @@
end subroutine mpas_tracer_advection_coefficients!}}}
- subroutine mpas_tracer_advection_tend(tend, s_old, s_new, diag, grid, dt, dminfo, cellsToSend, cellsToRecv)!{{{
+! 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)!{{{
- 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(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
type (mesh_type), intent(in) :: grid
- real (kind=RKIND) :: dt
+! real (kind=RKIND) :: dt
- type (dm_info), intent(in) :: dminfo
- type (exchange_list), pointer :: cellsToSend, cellsToRecv
+! 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
+! if(monotonicOn) then
+! call mpas_tracer_advection_mono_tend( tend, s_old, s_new, diag, grid, dt, dminfo, cellsToSend, cellsToRecv)
+! else
+ call mpas_tracer_advection_std_tend(tracers, uh, w, zDistance, zWeightUp, zWeightDown, grid, tend)
+! endif
end subroutine!}}}
subroutine mpas_tracer_advection_init(err)!{{{
Added: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_helpers.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_helpers.F (rev 0)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_helpers.F 2012-01-18 22:43:25 UTC (rev 1389)
@@ -0,0 +1,21 @@
+module mpas_tracer_advection_helpers
+
+ use mpas_kind_types
+
+ implicit none
+ save
+
+ contains
+
+ real function flux4(q_im2, q_im1, q_i, q_ip1, u)!{{{
+ real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, u
+ flux4 = u*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
+! flux4(q_im2, q_im1, q_i, q_ip1, ua) = ua*( 7.*(q_i + q_im1) - (q_ip1 + q_im2) )/12.0
+ end function!}}}
+
+ real function flux3( q_im2, q_im1, q_i, q_ip1, u, coef)!{{{
+ real (kind=RKIND) :: q_im2, q_im1, q_i, q_ip1, u, coef
+ flux3 = (u * (7.0 * (q_i + q_im1) - (q_ip1 + q_im2)) + coef * abs(u) * ((q_ip1 - q_im2) - 3.0*(q_i-q_im1)))/12.0
+ end function!}}}
+
+end module mpas_tracer_advection_helpers
Modified: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std.F 2012-01-18 20:39:35 UTC (rev 1388)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std.F 2012-01-18 22:43:25 UTC (rev 1389)
@@ -1,16 +1,23 @@
module mpas_tracer_advection_std
+ use mpas_kind_types
use mpas_grid_types
use mpas_configure
use mpas_dmpar
+ use mpas_tracer_advection_std_hadv
+ use mpas_tracer_advection_std_vadv
+
implicit none
private
save
+ public :: mpas_tracer_advection_std_tend, &
+ mpas_tracer_advection_std_init
+
contains
- subroutine mpas_tracer_advection_std_tend( tend, s_old, s_new, diag, grid, dt)!{{{
+ subroutine mpas_tracer_advection_std_tend(tracers, uh, w, zDistance, zWeightUp, zWeightDown, grid, tend)!{{{
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input: s - current model state
@@ -19,364 +26,29 @@
! Output: tend - computed scalar tendencies
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- 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(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
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
+ call mpas_tracer_advection_std_hadv_tend(tracers, uh, grid, tend)
+ call mpas_tracer_advection_std_vadv_tend(tracers, w, zDistance, zWeightUp, zWeightDown, grid, tend)
- 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 )!{{{
+ subroutine mpas_tracer_advection_std_init(err)!{{{
+ integer, intent(inout) :: err
- implicit none
- type (mesh_type) :: grid
+ integer :: err_tmp
- 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
+ err = 0
- integer :: cell1, cell2, iEdge, n, i, j, j_in, iCell
- integer :: cell_list(20), ordered_cell_list(20)
- logical :: addcell
+ call mpas_tracer_advection_std_hadv_init(err_tmp)
+ err = ior(err, err_tmp)
+ call mpas_tracer_advection_std_vadv_init(err_tmp)
+ err = ior(err, err_tmp)
- 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
+ end subroutine mpas_tracer_advection_std_init!}}}
- 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
Added: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_hadv.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_hadv.F (rev 0)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_hadv.F 2012-01-18 22:43:25 UTC (rev 1389)
@@ -0,0 +1,105 @@
+module mpas_tracer_advection_std_hadv
+
+ use mpas_kind_types
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ use mpas_tracer_advection_helpers
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_std_hadv_tend, &
+ mpas_tracer_advection_std_hadv_init
+
+ real (kind=RKIND) :: coef_3rd_order
+
+ logical :: hadvOn
+
+ contains
+
+ subroutine mpas_tracer_advection_std_hadv_tend(tracers, uh, grid, tend)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tracer tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: uh
+ type (mesh_type), intent(in) :: grid
+
+ integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
+ real (kind=RKIND) :: flux, tracer_weight
+
+ real (kind=RKIND), dimension(:), pointer :: 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( s_old % num_tracers, grid % nVertLevels ) :: flux_arr
+ real (kind=RKIND), dimension(:,:), allocatable :: flux_arr
+ integer :: nVertLevels, num_tracers
+
+ if (.not. hadvOn) return
+
+ cellsOnEdge => grid % cellsOnEdge % array
+ areaCell => grid % areaCell % 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
+ num_tracers = size(tracers, dim=1)
+
+ allocate(flux_arr(num_tracers, grid % nVertLevels))
+
+ ! horizontal flux divergence, accumulate in tracer_tend
+
+ 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
+ flux_arr(:,:) = 0.
+ do i=1,nAdvCellsForEdge(iEdge)
+ iCell = advCellsForEdge(i,iEdge)
+ do k=1,grid % nVertLevels
+ tracer_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uh(k,iEdge))*adv_coefs_3rd(i,iEdge)
+ do iTracer=1,num_tracers
+ flux_arr(iTracer,k) = flux_arr(iTracer,k) + tracer_weight* tracers(iTracer,k,iCell)
+ end do
+ end do
+ end do
+
+ do k=1,grid % nVertLevels
+ do iTracer=1,num_tracers
+ tend(iTracer,k,cell1) = tend(iTracer,k,cell1) - uh(k,iEdge)*flux_arr(iTracer,k)/areaCell(cell1)
+ tend(iTracer,k,cell2) = tend(iTracer,k,cell2) + uh(k,iEdge)*flux_arr(iTracer,k)/areaCell(cell2)
+ end do
+ end do
+ end if
+ end do
+
+ deallocate(flux_arr)
+
+ end subroutine mpas_tracer_advection_std_hadv_tend!}}}
+
+ subroutine mpas_tracer_advection_std_hadv_init(err)!{{{
+ integer, intent(inout) :: err
+
+ err = 0
+
+ hadvOn =.true.
+
+ coef_3rd_order = config_coef_3rd_order
+ end subroutine mpas_tracer_advection_std_hadv_init!}}}
+
+end module mpas_tracer_advection_std_hadv
Added: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_vadv.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_vadv.F (rev 0)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_vadv.F 2012-01-18 22:43:25 UTC (rev 1389)
@@ -0,0 +1,66 @@
+module mpas_tracer_advection_std_vadv
+
+ use mpas_kind_types
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ use mpas_tracer_advection_std_vadv2
+ use mpas_tracer_advection_std_vadv3
+ use mpas_tracer_advection_std_vadv4
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_std_vadv_tend, &
+ mpas_tracer_advection_std_vadv_init
+
+ logical :: order2, order3, order4
+
+ contains
+
+ subroutine mpas_tracer_advection_std_vadv_tend(tracers, w, zDistance, zWeightUp, zWeightDown, grid, tend)!{{{
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: w
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightUp, zWeightDown
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND) :: dt
+
+ if(order2) then
+ call mpas_tracer_advection_std_vadv2_tend(tracers, w, zDistance, zWeightUp, zWeightDown, grid, tend)
+ else if(order3) then
+ call mpas_tracer_advection_std_vadv3_tend(tracers, w, zDistance, zWeightUp, zWeightDown, grid, tend)
+ else if(order4) then
+ call mpas_tracer_advection_std_vadv4_tend(tracers, w, zDistance, zWeightUp, zWeightDown, grid, tend)
+ endif
+
+ end subroutine mpas_tracer_advection_std_vadv_tend!}}}
+
+ subroutine mpas_tracer_advection_std_vadv_init(err)!{{{
+ integer, intent(inout) :: err
+
+ err = 0
+
+ order2 = .false.
+ order3 = .false.
+ order4 = .false.
+
+ if (config_tracer_vadv_order == 2) then
+ order2 = .true.
+ else if (config_tracer_vadv_order == 3) then
+ order3 = .true.
+ else if (config_tracer_vadv_order == 4) then
+ order4 = .true.
+ else
+ print *, 'invalid value for config_tracer_vadv_order'
+ print *, 'options are 2, 3, or 4'
+ err = 1
+ endif
+
+ end subroutine mpas_tracer_advection_std_vadv_init!}}}
+
+end module mpas_tracer_advection_std_vadv
+
Added: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_vadv2.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_vadv2.F (rev 0)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_vadv2.F 2012-01-18 22:43:25 UTC (rev 1389)
@@ -0,0 +1,77 @@
+module mpas_tracer_advection_std_vadv2
+
+ use mpas_kind_types
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ use mpas_tracer_advection_helpers
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_std_vadv2_tend
+
+ contains
+
+ subroutine mpas_tracer_advection_std_vadv2_tend(tracers, w, zDistance, zWeightUp, zWeightDown, grid, tend)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tracer tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: w
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightUp, zWeightDown
+ type (mesh_type), intent(in) :: grid
+
+ integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
+ real (kind=RKIND) :: flux, tracer_edge, tracer_weight
+ real (kind=RKIND) :: tracer_weight_cell1, tracer_weight_cell2
+
+
+! real (kind=RKIND), dimension( s_old % num_tracers, grid % nVertLevels + 1 ) :: vert_flux
+ real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
+ integer :: nVertLevels, num_tracers
+
+ integer, parameter :: hadv_opt = 2
+
+ nVertLevels = grid % nVertLevels
+ num_tracers = size(tracers, dim=1)
+
+ allocate(vert_flux(num_tracers, nVertLevels+1))
+
+ !
+ ! vertical flux divergence
+ !
+
+ ! zero fluxes at top and bottom
+
+ vert_flux(:,1) = 0.
+ vert_flux(:,grid % nVertLevels+1) = 0.
+
+ do iCell=1,grid % nCellsSolve
+ do k = 2, nVertLevels
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = w(k,iCell)*(zWeightDown(k)*tracers(iTracer,k,iCell)+zWeightUp(k)*tracers(iTracer,k-1,iCell))
+ end do
+ end do
+ do k=1,grid % nVertLevelsSolve
+ do iTracer=1,num_tracers
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) - ( vert_flux(iTracer, k+1) - vert_flux(iTracer, k)) / zDistance(k)
+! tracers(iTracer,k,iCell) = ( tracers(iTracer,k,iCell)*h_old(k,iCell) &
+! + dt*( tracer_tend(iTracer,k,iCell) -rdnw(k)*(vert_flux(iTracer,k+1)-vert_flux(iTracer,k)) ) )/h_new(k,iCell)
+ end do
+ end do
+ end do
+
+ deallocate(vert_flux)
+
+ end subroutine mpas_tracer_advection_std_vadv2_tend!}}}
+
+end module mpas_tracer_advection_std_vadv2
Added: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_vadv3.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_vadv3.F (rev 0)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_vadv3.F 2012-01-18 22:43:25 UTC (rev 1389)
@@ -0,0 +1,89 @@
+module mpas_tracer_advection_std_vadv3
+
+ use mpas_kind_types
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ use mpas_tracer_advection_helpers
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_std_vadv3_tend
+
+ real (kind=RKIND) :: coef_3rd_order
+
+ contains
+
+ subroutine mpas_tracer_advection_std_vadv3_tend(tracers, w, zDistance, zWeightUp, zWeightDown, grid, tend)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tracer tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: w
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightUp, zWeightDown
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND) :: dt
+
+ integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
+
+! real (kind=RKIND), dimension( s_old % num_tracers, grid % nVertLevels + 1 ) :: vert_flux
+ real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
+ integer :: nVertLevels, num_tracers
+
+ integer, parameter :: hadv_opt = 2
+
+ coef_3rd_order = config_coef_3rd_order
+
+ nVertLevels = grid % nVertLevels
+ num_tracers = size(tracers, dim=1)
+
+ allocate(vert_flux(num_tracers, nVertLevels+1))
+
+ vert_flux(:,1) = 0.
+ vert_flux(:,grid % nVertLevels+1) = 0.
+
+ do iCell=1,grid % nCellsSolve
+
+ k = 2
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = w(k,iCell)*(zWeightDown(k)*tracers(iTracer,k,iCell)+zWeightUp(k)*tracers(iTracer,k-1,iCell))
+ enddo
+
+ do k=3,nVertLevels-1
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = flux3( tracers(iTracer,k-2,iCell),tracers(iTracer,k-1,iCell), &
+ tracers(iTracer,k ,iCell),tracers(iTracer,k+1,iCell), &
+ w(k,iCell), coef_3rd_order )
+ end do
+ end do
+
+ k = nVertLevels
+
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = w(k,iCell)*(zWeightDown(k)*tracers(iTracer,k,iCell)+zWeightUp(k)*tracers(iTracer,k-1,iCell))
+ enddo
+
+ do k=1,grid % nVertLevelsSolve
+ do iTracer=1,num_tracers
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) - (vert_flux(iTracer, k+1) - vert_flux(iTracer, k)) / zDistance(k)
+
+! tracers(iTracer,k,iCell) = ( tracers(iTracer,k,iCell)*h_old(k,iCell) &
+! + dt*( tracer_tend(iTracer,k,iCell) -rdnw(k)*(vert_flux(iTracer,k+1)-vert_flux(iTracer,k)) ) )/h_new(k,iCell)
+ end do
+ end do
+ end do
+
+ deallocate(vert_flux)
+
+ end subroutine mpas_tracer_advection_std_vadv3_tend!}}}
+
+end module mpas_tracer_advection_std_vadv3
Added: branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_vadv4.F
===================================================================
--- branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_vadv4.F (rev 0)
+++ branches/ocean_projects/advection_routines/src/operators/mpas_tracer_advection_std_vadv4.F 2012-01-18 22:43:25 UTC (rev 1389)
@@ -0,0 +1,85 @@
+module mpas_tracer_advection_std_vadv4
+
+ use mpas_kind_types
+ use mpas_grid_types
+ use mpas_configure
+ use mpas_dmpar
+
+ use mpas_tracer_advection_helpers
+
+ implicit none
+ private
+ save
+
+ public :: mpas_tracer_advection_std_vadv4_tend
+
+ contains
+
+ subroutine mpas_tracer_advection_std_vadv4_tend(tracers, w, zDistance, zWeightUp, zWeightDown, grid, tend)!{{{
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !
+ ! Input: s - current model state
+ ! grid - grid metadata
+ !
+ ! Output: tend - computed tracer tendencies
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
+ real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+ real (kind=RKIND), dimension(:,:), intent(in) :: w
+ real (kind=RKIND), dimension(:), intent(in) :: zDistance, zWeightUp, zWeightDown
+ type (mesh_type), intent(in) :: grid
+ real (kind=RKIND) :: dt
+
+ integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
+
+! real (kind=RKIND), dimension( s_old % num_tracers, grid % nVertLevels + 1 ) :: vert_flux
+ real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
+ integer :: nVertLevels, num_tracers
+
+ nVertLevels = grid % nVertLevels
+ num_tracers = size(tracers, dim=1)
+
+ allocate(vert_flux(num_tracers, nVertLevels+1))
+
+ ! vertical flux divergence
+ !
+
+ ! zero fluxes at top and bottom
+
+ vert_flux(:,1) = 0.
+ vert_flux(:,grid % nVertLevels+1) = 0.
+
+ do iCell=1,grid % nCellsSolve
+
+ k = 2
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = w(k,iCell)*(zWeightDown(k)*tracers(iTracer,k,iCell)+zWeightUp(k)*tracers(iTracer,k-1,iCell))
+ enddo
+ do k=3,nVertLevels-1
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = flux4( tracers(iTracer,k-2,iCell),tracers(iTracer,k-1,iCell), &
+ tracers(iTracer,k ,iCell),tracers(iTracer,k+1,iCell), w(k,iCell) )
+ end do
+ end do
+
+ k = nVertLevels
+ do iTracer=1,num_tracers
+ vert_flux(iTracer,k) = w(k,iCell)*(zWeightDown(k)*tracers(iTracer,k,iCell)+zWeightUp(k)*tracers(iTracer,k-1,iCell))
+ enddo
+
+ do k=1,grid % nVertLevelsSolve
+ do iTracer=1,num_tracers
+ tend(iTracer, k, iCell) = tend(iTracer, k, iCell) - (vert_flux(iTracer, k+1) - vert_flux(iTracer, k)) / zDistance(k)
+! tracers(iTracer,k,iCell) = ( tracers(iTracer,k,iCell)*h_old(k,iCell) &
+! + dt*( tracer_tend(iTracer,k,iCell) -rdnw(k)*(vert_flux(iTracer,k+1)-vert_flux(iTracer,k)) ) )/h_new(k,iCell)
+ end do
+ end do
+
+ end do
+
+ deallocate(vert_flux)
+
+ end subroutine mpas_tracer_advection_std_vadv4_tend!}}}
+
+end module mpas_tracer_advection_std_vadv4
</font>
</pre>