<p><b>dwj07@fsu.edu</b> 2012-03-09 13:42:24 -0700 (Fri, 09 Mar 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Cleaning up namelist/registry variables.<br>
Adding doxygen comments to mpas_ocn_tendency routines.<br>
<br>
Removing shared advection routines, since they will temporarily be placed in core_ocean.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/monotonic_advection/namelist.input.ocean
===================================================================
--- branches/ocean_projects/monotonic_advection/namelist.input.ocean 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/namelist.input.ocean 2012-03-09 20:42:24 UTC (rev 1614)
@@ -76,10 +76,8 @@
&advection
config_vert_tracer_adv = 'stencil'
config_vert_tracer_adv_order = 2
- config_tracer_adv_order = 2
- config_tracer_vadv_order = 2
+ config_horiz_tracer_adv_order = 2
config_thickness_adv_order = 2
- config_positive_definite = .false.
config_monotonic = .false.
/
&restore
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/Registry 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/Registry 2012-03-09 20:42:24 UTC (rev 1614)
@@ -66,12 +66,10 @@
namelist real vmix_tanh config_zWidth_tanh 100
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_vert_tracer_adv_order 4
+namelist integer advection config_horiz_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 advection config_check_monotonicity false
namelist logical restore config_restoreTS false
@@ -180,9 +178,6 @@
var persistent real referenceBottomDepth ( nVertLevels ) 0 iro referenceBottomDepth mesh - -
var persistent real referenceBottomDepthTopOfCell ( nVertLevelsP1 ) 0 - referenceBottomDepthTopOfCell mesh - -
var persistent real hZLevel ( nVertLevels ) 0 iro hZLevel mesh - -
-var persistent real hMeanTopZLevel ( nVertLevels ) 0 - hMeanTopZLevel mesh - -
-var persistent real hRatioZLevelK ( nVertLevels ) 0 - hRatioZLevelK mesh - -
-var persistent real hRatioZLevelKm1 ( nVertLevels ) 0 - hRatioZLevelKm1 mesh - -
% Boundary conditions: read from input, saved in restart and written to output
var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_mpas_core.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_mpas_core.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -278,29 +278,13 @@
:block % mesh % nVertLevels,iCell) = 0.0
! mrp changed to 0
! :block % mesh % nVertLevels,iCell) = -1e34
-
-! mrp 110516, added just to test for conservation of tracers
- block % state % time_levs(1) % state % tracers % array(block % state % time_levs(1) % state % index_tracer1,:,iCell) = 1.0
-
end do
- if (.not. config_do_restart) then
-
-! mrp 110808 add, so that variables are copied to * variables for split explicit
- do i=2,nTimeLevs
- call mpas_copy_state(block % state % time_levs(i) % state, &
+ do i=2,nTimeLevs
+ call mpas_copy_state(block % state % time_levs(i) % state, &
block % state % time_levs(1) % state)
- end do
-! mrp 110808 add end
+ end do
-
- else
- do i=2,nTimeLevs
- call mpas_copy_state(block % state % time_levs(i) % state, &
- block % state % time_levs(1) % state)
- end do
- endif
-
end subroutine mpas_init_block!}}}
subroutine mpas_core_run(domain, output_obj, output_frame)!{{{
@@ -514,7 +498,7 @@
integer :: iTracer, cell, cell1, cell2
real (kind=RKIND) :: uhSum, hSum, hEdge1
- real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth, referenceBottomDepthTopOfCell, hMeanTopZLevel, hRatioZLevelK, hRatioZLevelKm1
+ real (kind=RKIND), dimension(:), pointer :: referenceBottomDepth, referenceBottomDepthTopOfCell
real (kind=RKIND), dimension(:,:), pointer :: h
integer :: nVertLevels
@@ -527,9 +511,6 @@
referenceBottomDepth => block % mesh % referenceBottomDepth % array
referenceBottomDepthTopOfCell => block % mesh % referenceBottomDepthTopOfCell % array
nVertLevels = block % mesh % nVertLevels
- hMeanTopZLevel => block % mesh % hMeanTopZLevel % array
- hRatioZLevelK => block % mesh % hRatioZLevelK % array
- hRatioZLevelKm1 => block % mesh % hRatioZLevelKm1 % array
! mrp 120208 right now hZLevel is in the grid.nc file.
! We would like to transition to using referenceBottomDepth
@@ -537,14 +518,8 @@
! When the transition is done, hZLevel can be removed from
! registry and the following four lines deleted.
referenceBottomDepth(1) = block % mesh % hZLevel % array(1)
- hMeanTopZLevel(1) = 0.0
- hRatioZLevelK(1) = 0.0
- hRatioZLevelKm1(1) = 0.0
do k = 2,nVertLevels
referenceBottomDepth(k) = referenceBottomDepth(k-1) + block % mesh % hZLevel % array(k)
- hMeanTopZLevel(k) = 0.5*(block % mesh % hZLevel % array(k-1) + block % mesh % hZLevel % array(k))
- hRatioZLevelK(k) = 0.5*block % mesh % hZLevel % array(k)/hMeanTopZLevel(k)
- hRatioZLevelKm1(k) = 0.5*block % mesh % hZLevel % array(k-1)/hMeanTopZLevel(k)
end do
! TopOfCell needed where zero depth for the very top may be referenced.
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tendency.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tendency.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -99,22 +99,12 @@
!
!-----------------------------------------------------------------------
- subroutine ocn_tend_h(tend, s, d, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute height and normal wind tendencies, as well as diagnostic variables
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed tendencies for prognostic variables
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+ subroutine ocn_tend_h(tend, s, grid)!{{{
implicit none
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (diagnostics_type), intent(in) :: d
- type (mesh_type), intent(in) :: grid
+ type (tend_type), intent(inout) :: tend !< Input/Output: Tendency structure
+ type (state_type), intent(in) :: s !< Input: State information
+ type (mesh_type), intent(in) :: grid !< Input: Grid information
real (kind=RKIND), dimension(:,:), pointer :: h_edge, u, wTop, tend_h
@@ -168,21 +158,12 @@
!-----------------------------------------------------------------------
subroutine ocn_tend_u(tend, s, d, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute height and normal wind tendencies, as well as diagnostic variables
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed tendencies for prognostic variables
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
implicit none
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (diagnostics_type), intent(in) :: d
- type (mesh_type), intent(in) :: grid
+ type (tend_type), intent(inout) :: tend !< Input/Output: Tendency structure
+ type (state_type), intent(in) :: s !< Input: State information
+ type (diagnostics_type), intent(in) :: d !< Input: Diagnostic information
+ type (mesh_type), intent(in) :: grid !< Input: Grid information
real (kind=RKIND), dimension(:,:), pointer :: &
h_edge, h, u, rho, zMid, pressure, &
@@ -288,27 +269,15 @@
!> This routine computes the scalar (tracer) tendency for the ocean
!
!-----------------------------------------------------------------------
-
subroutine ocn_tend_scalar(tend, s, d, grid, dt)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! Input: s - current model state
- ! grid - grid metadata
- ! note: the variable s % tracers really contains the tracers,
- ! not tracers*h
- !
- ! Output: tend - computed scalar tendencies
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
implicit none
- type (tend_type), intent(inout) :: tend
- type (state_type), intent(in) :: s
- type (diagnostics_type), intent(in) :: d
- type (mesh_type), intent(in) :: grid
- real (kind=RKIND), intent(in) :: dt
+ type (tend_type), intent(inout) :: tend !< Input/Output: Tendency structure
+ type (state_type), intent(in) :: s !< Input: State information
+ type (diagnostics_type), intent(in) :: d !< Input: Diagnostic information
+ type (mesh_type), intent(in) :: grid !< Input: Grid information
+ real (kind=RKIND), intent(in) :: dt !< Input: Time step
- real (kind=RKIND), dimension(:), pointer :: hRatioZLevelK, hRatioZLevelKm1, zTopZLevel, hMeanTopZLevel
real (kind=RKIND), dimension(:,:), pointer :: &
u, h,wTop, h_edge, vertDiffTopOfCell, tend_h, uh
real (kind=RKIND), dimension(:,:,:), pointer :: &
@@ -325,10 +294,6 @@
h_edge => s % h_edge % array
vertDiffTopOfCell => d % vertDiffTopOfCell % array
- hRatioZLevelK => grid % hRatioZLevelK % array
- hRatioZLevelKm1 => grid % hRatioZLevelKm1 % array
- hMeanTopZLevel => grid % hMeanTopZLevel % array
-
tend_tr => tend % tracers % array
tend_h => tend % h % array
@@ -418,19 +383,11 @@
!-----------------------------------------------------------------------
subroutine ocn_diagnostic_solve(dt, s, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute diagnostic fields used in the tendency computations
- !
- ! Input: grid - grid metadata
- !
- ! Output: s - computed diagnostics
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
implicit none
- real (kind=RKIND), intent(in) :: dt
- type (state_type), intent(inout) :: s
- type (mesh_type), intent(in) :: grid
+ real (kind=RKIND), intent(in) :: dt !< Input: Time step
+ type (state_type), intent(inout) :: s !< Input/Output: State information
+ type (mesh_type), intent(in) :: grid !< Input: Grid information
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j
@@ -887,23 +844,13 @@
!> This routine computes the vertical velocity in the top layer for the ocean
!
!-----------------------------------------------------------------------
-
subroutine ocn_wtop(s1,s2, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute diagnostic fields used in the tendency computations
- !
- ! Input: grid - grid metadata
- !
- ! Output: s - computed diagnostics
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
implicit none
- type (state_type), intent(inout) :: s1
- type (state_type), intent(inout) :: s2
- type (mesh_type), intent(in) :: grid
+ type (state_type), intent(inout) :: s1 !< Input/Output: State 1 information
+ type (state_type), intent(inout) :: s2 !< Input/Output: State 2 information
+ type (mesh_type), intent(in) :: grid !< Input: Grid information
- ! mrp 110512 could clean this out, remove pointers?
integer :: iEdge, iCell, iVertex, k, cell1, cell2, vertex1, vertex2, eoe, i, j, cov
real (kind=RKIND) :: flux, vorticity_abs, h_vertex, workpv, rho0Inv, hSum
@@ -1077,19 +1024,10 @@
!-----------------------------------------------------------------------
subroutine ocn_fuperp(s, grid)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Put f*uBcl^{perp} in u as a work variable
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed tendencies for prognostic variables
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
implicit none
- type (state_type), intent(inout) :: s
- type (mesh_type), intent(in) :: grid
+ type (state_type), intent(inout) :: s !< Input/Output: State information
+ type (mesh_type), intent(in) :: grid !< Input: Grid information
! mrp 110512 I just split compute_tend into compute_tend_u and compute_tend_h.
! Some of these variables can be removed, but at a later time.
@@ -1152,9 +1090,8 @@
!> other tendency routines.
!
!-----------------------------------------------------------------------
-
subroutine ocn_tendency_init(err)!{{{
- integer, intent(out) :: err
+ integer, intent(out) :: err !< Output: Error flag
err = 0
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -173,7 +173,7 @@
if (.not.config_implicit_vertical_mix) then
call ocn_vmix_coefs(block % mesh, provis, block % diagnostics, err)
end if
- call ocn_tend_h(block % tend, provis, block % diagnostics, block % mesh)
+ call ocn_tend_h(block % tend, provis, block % mesh)
call ocn_tend_u(block % tend, provis, block % diagnostics, block % mesh)
! mrp 110718 filter btr mode out of u_tend
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_split.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_split.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -200,7 +200,7 @@
if (.not.config_implicit_vertical_mix) then
call ocn_vmix_coefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
end if
- call ocn_tend_u(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+ call ocn_tend_u(block % tend, block % state % time_levs(2) % state, block % diagnostics, block % mesh)
block => block % next
end do
@@ -721,7 +721,7 @@
do while (associated(block))
call ocn_wtop(block % state % time_levs(1) % state,block % state % time_levs(2) % state, block % mesh)
- call ocn_tend_h (block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh)
+ call ocn_tend_h(block % tend, block % state % time_levs(2) % state, block % mesh)
block => block % next
end do
@@ -738,7 +738,7 @@
block => domain % blocklist
do while (associated(block))
- call ocn_tend_scalar(block % tend, block % state % time_levs(2) % state , block % diagnostics, block % mesh, dt)
+ call ocn_tend_scalar(block % tend, block % state % time_levs(2) % state, block % diagnostics, block % mesh, dt)
block => block % next
end do
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_std_vadv.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_std_vadv.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_std_vadv.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -85,11 +85,11 @@
order3 = .false.
order4 = .false.
- if (config_tracer_vadv_order == 2) then
+ if (config_vert_tracer_adv_order == 2) then
order2 = .true.
- else if (config_tracer_vadv_order == 3) then
+ else if (config_vert_tracer_adv_order == 3) then
order3 = .true.
- else if (config_tracer_vadv_order == 4) then
+ else if (config_vert_tracer_adv_order == 4) then
order4 = .true.
else
print *, 'invalid value for config_tracer_vadv_order'
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_hadv2.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_hadv2.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_hadv2.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -182,7 +182,7 @@
err = 0
hadv2On = .false.
- if (config_tracer_adv_order == 2) then
+ if (config_horiz_tracer_adv_order == 2) then
hadv2On = .true.
end if
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_hadv3.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_hadv3.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_hadv3.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -221,11 +221,10 @@
err = 0
hadv3On = .false.
- if (config_tracer_adv_order == 3) then
+ if (config_horiz_tracer_adv_order == 3) then
hadv3On = .true.
- coef_3rd_order = 1.0
- if (config_monotonic) coef_3rd_order = 0.25
+ coef_3rd_order = config_coef_3rd_order
end if
!--------------------------------------------------------------------
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_hadv4.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_hadv4.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_hadv4.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -213,7 +213,7 @@
err = 0
hadv4On = .false.
- if (config_tracer_adv_order == 4) then
+ if (config_horiz_tracer_adv_order == 4) then
hadv4On = .true.
end if
Modified: branches/ocean_projects/monotonic_advection/src/operators/Makefile
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/Makefile 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/operators/Makefile 2012-03-09 20:42:24 UTC (rev 1614)
@@ -2,16 +2,7 @@
OBJS = mpas_rbf_interpolation.o \
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 \
- mpas_tracer_advection_std_vadv2.o \
- mpas_tracer_advection_std_vadv3.o \
- mpas_tracer_advection_std_vadv4.o \
- mpas_tracer_advection_helpers.o
+ mpas_spline_interpolation.o
all: operators
@@ -24,24 +15,6 @@
mpas_spline_interpolation:
-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
-
-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
Deleted: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -1,217 +0,0 @@
-module mpas_tracer_advection
-
- use mpas_kind_types
- use mpas_grid_types
- use mpas_configure
-
- 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, dimension(:), pointer :: cell_list, ordered_cell_list
- integer :: cell1, cell2, iEdge, n, i, j, j_in, iCell
- 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
-
- allocate(cell_list(grid % maxEdges2 + 2))
- allocate(ordered_cell_list(grid % maxEdges2 + 2))
-
- 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
-
- deallocate(cell_list)
- deallocate(ordered_cell_list)
-
- end subroutine mpas_tracer_advection_coefficients!}}}
-
- subroutine mpas_tracer_advection_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend)!{{{
-
- real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
- real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
- real (kind=RKIND), dimension(:,:), intent(in) :: uh, w, h
- real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !< Input: Distance between vertical interfaces of a cell
- real (kind=RKIND), intent(in) :: dt
- type (mesh_type), intent(in) :: grid
- real (kind=RKIND), dimension(:,:), intent(in) :: tend_h
-
-
- if(monotonicOn) then
- call mpas_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend)
- else
- call mpas_tracer_advection_std_tend(tracers, uh, w, verticalCellSize, grid, tend)
- endif
- end subroutine mpas_tracer_advection_tend!}}}
-
- subroutine mpas_tracer_advection_init(err)!{{{
-
- integer, intent(inout) :: err
-
- integer :: err_tmp
-
- err = 0
-
- call mpas_tracer_advection_std_init(err_tmp)
-
- err = ior(err, err_tmp)
-
- monotonicOn = .false.
-
- if(config_monotonic) then
- monotonicOn = .true.
- endif
-
- end subroutine mpas_tracer_advection_init!}}}
-
-end module mpas_tracer_advection
Deleted: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_helpers.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_helpers.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_helpers.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -1,23 +0,0 @@
-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), intent(in) :: q_im2, q_im1, q_i, q_ip1, u
- flux4 = u*( 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), intent(in) :: q_im2, q_im1, q_i, q_ip1, u, coef
-
- !dwj 02/21/12 flux3 is different in ocean and atmosphere
- !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
- 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
Deleted: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_mono.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_mono.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -1,323 +0,0 @@
-module mpas_tracer_advection_mono
-
- 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_mono_tend
-
- contains
-
- subroutine mpas_tracer_advection_mono_tend(tracers, uh, w, h, verticalCellSize, dt, grid, tend_h, tend)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers !< Input: current tracer values
- real (kind=RKIND), dimension(:,:), intent(in) :: uh !< Input: Thichness weighted velocitiy
- real (kind=RKIND), dimension(:,:), intent(in) :: w !< Input: Vertical velocitiy
- real (kind=RKIND), dimension(:,:), intent(in) :: h !< Input: Thickness
- real (kind=RKIND), dimension(:,:), intent(in) :: verticalCellSize !< Input: Distance between vertical interfaces of a cell
- real (kind=RKIND), dimension(:,:), intent(in) :: tend_h !< Input: Tendency for thickness field
- real (kind=RKIND), intent(in) :: dt !< Input: Timestep
- type (mesh_type), intent(in) :: grid
- real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
-
- integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
- integer :: nVertLevels, num_tracers, nCells, nEdges, nCellsSolve
- integer, dimension(:), pointer :: nAdvCellsForEdge, maxLevelCell, maxLevelEdgeTop, nEdgesOnCell
- integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnCell, advCellsForEdge
-
- real (kind=RKIND) :: coef_3rd_order, flux_upwind, tracer_min_new, tracer_max_new, tracer_upwind_new, scale_factor
- real (kind=RKIND) :: flux, tracer_weight, invDvEdge, invAreaCell1, invAreaCell2
- real (kind=RKIND) :: cur_max, cur_min, new_max, new_min
- real (kind=RKIND) :: verticalWeightK, verticalWeightKm1
- real (kind=RKIND), dimension(:), pointer :: dvEdge, areaCell
- real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_3rd
- real (kind=RKIND), dimension(:,:), pointer :: tracer_cur, tracer_new, upwind_tendency, inv_h_new, tracer_max, tracer_min
- real (kind=RKIND), dimension(:,:), pointer :: flux_incoming, flux_outgoing, high_order_horiz_flux, high_order_vert_flux
-
- real (kind=RKIND), parameter :: eps = 1.e-10
-
- coef_3rd_order = config_coef_3rd_order
-
- ! Initialize pointers
- dvEdge => grid % dvEdge % array
- cellsOnEdge => grid % cellsOnEdge % array
- cellsOnCell => grid % cellsOnCell % array
- areaCell => grid % areaCell % array
-
- nEdgesOnCell => grid % nEdgesOnCell % 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
-
- nCells = grid % nCells
- nCellsSolve = grid % nCellsSolve
- nEdges = grid % nEdges
- nVertLevels = grid % nVertLevels
- num_tracers = size(tracers,dim=1)
-
- ! allocate nCells arrays
-
- allocate(tracer_new(nVertLevels, nCells))
- allocate(tracer_cur(nVertLevels, nCells))
- allocate(upwind_tendency(nVertLevels, nCells))
- allocate(inv_h_new(nVertLevels, nCells))
- allocate(tracer_max(nVertLevels, nCells))
- allocate(tracer_min(nVertLevels, nCells))
- allocate(flux_incoming(nVertLevels, nCells))
- allocate(flux_outgoing(nVertLevels, nCells))
-
- ! allocate nEdges arrays
- allocate(high_order_horiz_flux(nVertLevels, nEdges))
-
- ! allocate nVertLevels+1 and nCells arrays
- allocate(high_order_vert_flux(nVertLevels+1, nCells))
-
- do iCell = 1, nCells
- do k=1, maxLevelCell(iCell)
- inv_h_new(k, iCell) = 1.0 / (h(k, iCell) + dt * tend_h(k, iCell))
- end do
- end do
-
- ! Loop over tracers. One tracer is advected at a time. It is copied into a temporary array in order to improve locality
- do iTracer = 1, num_tracers
- ! Initialize variables for use in this iTracer iteration
- do iCell = 1, nCells
- do k=1, maxLevelCell(iCell)
- tracer_cur(k,iCell) = tracers(iTracer,k,iCell)
- upwind_tendency(k, iCell) = 0.0
-
- !tracer_new is supposed to be the "new" tracer state. This allows bounds checks.
- if (config_check_monotonicity) then
- tracer_new(k,iCell) = 0.0
- end if
- end do ! k loop
- end do ! iCell loop
-
- high_order_vert_flux = 0.0
- high_order_horiz_flux = 0.0
-
- ! Compute the high order vertical flux. Also determine bounds on tracer_cur.
- do iCell = 1, nCells
- k = 1
- tracer_max(k,iCell) = max(tracer_cur(k,iCell),tracer_cur(k+1,iCell))
- tracer_min(k,iCell) = min(tracer_cur(k,iCell),tracer_cur(k+1,iCell))
-
- k = 2
- verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
- verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
- high_order_vert_flux(k,iCell) = w(k,iCell)*(verticalWeightK*tracer_cur(k,iCell)+verticalWeightKm1*tracer_cur(k-1,iCell))
- tracer_max(k,iCell) = max(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
- tracer_min(k,iCell) = min(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
-
- do k=3,maxLevelCell(iCell)-1
- high_order_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 )
- tracer_max(k,iCell) = max(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
- tracer_min(k,iCell) = min(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
- end do
-
- k = maxLevelCell(iCell)
- verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
- verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) + verticalCellSize(k-1, iCell))
- high_order_vert_flux(k,iCell) = w(k,iCell)*(verticalWeightK*tracer_cur(k,iCell)+verticalWeightKm1*tracer_cur(k-1,iCell))
- tracer_max(k,iCell) = max(tracer_cur(k,iCell),tracer_cur(k-1,iCell))
- tracer_min(k,iCell) = min(tracer_cur(k,iCell),tracer_cur(k-1,iCell))
-
- ! pull tracer_min and tracer_max from the (horizontal) surrounding cells
- do i = 1, nEdgesOnCell(iCell)
- do k=1, maxLevelCell(cellsOnCell(i, iCell))
- tracer_max(k,iCell) = max(tracer_max(k,iCell),tracer_cur(k, cellsOnCell(i,iCell)))
- tracer_min(k,iCell) = min(tracer_min(k,iCell),tracer_cur(k, cellsOnCell(i,iCell)))
- end do ! k loop
- end do ! i loop over nEdgesOnCell
- end do ! iCell Loop
-
- ! Compute the high order horizontal flux
- do iEdge = 1, nEdges
- do i = 1, nAdvCellsForEdge(iEdge)
- iCell = advCellsForEdge(i,iEdge)
- do k = 1, maxLevelCell(iCell)
- tracer_weight = uh(k,iEdge)*(adv_coefs(i,iEdge) + coef_3rd_order*sign(1.,uh(k,iEdge))*adv_coefs_3rd(i,iEdge))
- high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) + tracer_weight* tracer_cur(k,iCell)
- end do ! k loop
- end do ! i loop over nAdvCellsForEdge
- end do ! iEdge loop
-
- ! low order upwind vertical flux (monotonic and diffused)
- ! Remove low order flux from the high order flux.
- ! Store left over high order flux in high_order_vert_flux array.
- ! Upwind fluxes are accumulated in upwind_tendency
- do iCell = 1, nCells
- do k = 2, maxLevelCell(iCell)
- ! dwj 02/03/12 Ocean and Atmosphere are different in vertical
-! flux_upwind = max(0.,w(k,iCell))*tracer_cur(k-1,iCell) + min(0.,w(k,iCell))*tracer_cur(k,iCell)
- flux_upwind = min(0.,w(k,iCell))*tracer_cur(k-1,iCell) + max(0.,w(k,iCell))*tracer_cur(k,iCell)
- upwind_tendency(k-1,iCell) = upwind_tendency(k-1,iCell) + flux_upwind
- upwind_tendency(k ,iCell) = upwind_tendency(k ,iCell) - flux_upwind
- high_order_vert_flux(k,iCell) = high_order_vert_flux(k,iCell) - flux_upwind
- end do ! k loop
-
- ! flux_incoming contains the total remaining high order flux into iCell
- ! it is positive.
- ! flux_outgoing contains the total remaining high order flux out of iCell
- ! it is negative
- do k = 1, maxLevelCell(iCell)
- ! dwj 02/03/12 Ocean and Atmosphere are different in vertical
-! flux_incoming (k,iCell) = -(min(0.,high_order_vert_flux(k+1,iCell))-max(0.,high_order_vert_flux(k,iCell)))
-! flux_outgoing(k,iCell) = -(max(0.,high_order_vert_flux(k+1,iCell))-min(0.,high_order_vert_flux(k,iCell)))
-
- flux_incoming (k, iCell) = max(0.0, high_order_vert_flux(k+1, iCell)) - min(0.0, high_order_vert_flux(k, iCell))
- flux_outgoing(k, iCell) = min(0.0, high_order_vert_flux(k+1, iCell)) - max(0.0, high_order_vert_flux(k, iCell))
- end do ! k Loop
- end do ! iCell Loop
-
- ! low order upwind horizontal flux (monotinc and diffused)
- ! Remove low order flux from the high order flux
- ! Store left over high order flux in high_order_horiz_flux array
- ! Upwind fluxes are accumulated in upwind_tendency
- do iEdge = 1, nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
-
- 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))
- high_order_horiz_flux(k,iEdge) = high_order_horiz_flux(k,iEdge) - flux_upwind
-
- upwind_tendency(k,cell1) = upwind_tendency(k,cell1) - flux_upwind * invAreaCell1
- upwind_tendency(k,cell2) = upwind_tendency(k,cell2) + flux_upwind * invAreaCell2
-
- ! Accumulate remaining high order fluxes
- flux_outgoing(k,cell1) = flux_outgoing(k,cell1) - max(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell1
- flux_incoming (k,cell1) = flux_incoming (k,cell1) - min(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell1
- flux_outgoing(k,cell2) = flux_outgoing(k,cell2) + min(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell2
- flux_incoming (k,cell2) = flux_incoming (k,cell2) + max(0.,high_order_horiz_flux(k,iEdge)) * invAreaCell2
- end do ! k loop
- end do ! iEdge loop
-
- ! Build the factors for the FCT
- ! Computed using the bounds that were computed previously, and the bounds on the newly updated value
- ! Factors are placed in the flux_incoming and flux_outgoing arrays
- do iCell = 1, nCells
- do k = 1, maxLevelCell(iCell)
- tracer_min_new = (tracer_cur(k,iCell)*h(k,iCell) + dt*(upwind_tendency(k,iCell)+flux_outgoing(k,iCell))) * inv_h_new(k,iCell)
- tracer_max_new = (tracer_cur(k,iCell)*h(k,iCell) + dt*(upwind_tendency(k,iCell)+flux_incoming(k,iCell))) * inv_h_new(k,iCell)
- tracer_upwind_new = (tracer_cur(k,iCell)*h(k,iCell) + dt*upwind_tendency(k,iCell)) * inv_h_new(k,iCell)
-
- scale_factor = (tracer_max(k,iCell)-tracer_upwind_new)/(tracer_max_new-tracer_upwind_new+eps)
- flux_incoming(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
-
- scale_factor = (tracer_upwind_new-tracer_min(k,iCell))/(tracer_upwind_new-tracer_min_new+eps)
- flux_outgoing(k,iCell) = min( 1.0, max( 0.0, scale_factor) )
- end do ! k loop
- end do ! iCell loop
-
- ! rescale the high order horizontal fluxes
- do iEdge = 1, nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- do k = 1, maxLevelEdgeTop(iEdge)
- flux = high_order_horiz_flux(k,iEdge)
- flux = max(0.,flux) * min(flux_outgoing(k,cell1), flux_incoming(k,cell2)) &
- + min(0.,flux) * min(flux_incoming(k,cell1), flux_outgoing(k,cell2))
- high_order_horiz_flux(k,iEdge) = flux
- end do ! k loop
- end do ! iEdge loop
-
- ! rescale the high order vertical flux
- do iCell = 1, nCellsSolve
- do k = 2, maxLevelCell(iCell)
- flux = high_order_vert_flux(k,iCell)
- ! dwj 02/03/12 Ocean and Atmosphere are different in vertical.
-! flux = max(0.,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k ,iCell)) &
-! + min(0.,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell))
- flux = max(0.,flux) * min(flux_outgoing(k ,iCell), flux_incoming(k-1,iCell)) &
- + min(0.,flux) * min(flux_outgoing(k-1,iCell), flux_incoming(k ,iCell))
- high_order_vert_flux(k,iCell) = flux
- end do ! k loop
- end do ! iCell loop
-
- ! Accumulate the scaled high order horizontal tendencies
- do iEdge = 1, nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
-
- invAreaCell1 = 1.0 / areaCell(cell1)
- invAreaCell2 = 1.0 / areaCell(cell2)
- do k=1,maxLevelEdgeTop(iEdge)
- tend(iTracer, k, cell1) = tend(iTracer, k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
- tend(iTracer, k, cell2) = tend(iTracer, k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
-
- if (config_check_monotonicity) then
- !tracer_new holds a tendency for now.
- tracer_new(k, cell1) = tracer_new(k, cell1) - high_order_horiz_flux(k, iEdge) * invAreaCell1
- tracer_new(k, cell2) = tracer_new(k, cell2) + high_order_horiz_flux(k, iEdge) * invAreaCell2
- end if
- end do ! k loop
- end do ! iEdge loop
-
- ! Accumulate the scaled high order vertical tendencies, and the upwind tendencies
- do iCell = 1, nCellsSolve
- do k = 1,maxLevelCell(iCell)
- tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
-
- if (config_check_monotonicity) then
- !tracer_new holds a tendency for now. Only for a check on monotonicity
- tracer_new(k, iCell) = tracer_new(k, iCell) + (high_order_vert_flux(k+1, iCell) - high_order_vert_flux(k, iCell)) + upwind_tendency(k,iCell)
-
- !tracer_new is now the new state of the tracer. Only for a check on monotonicity
- tracer_new(k, iCell) = (tracer_cur(k, iCell)*h(k, iCell) + dt * tracer_new(k, iCell)) * inv_h_new(k, iCell)
- end if
- end do ! k loop
- end do ! iCell loop
-
- if (config_check_monotonicity) then
- !build min and max bounds on old and new tracer for check on monotonicity.
- cur_min = minval(tracer_cur(:,1:nCellsSolve))
- cur_max = maxval(tracer_cur(:,1:nCellsSolve))
- new_min = minval(tracer_new(:,1:nCellsSolve))
- new_max = maxval(tracer_new(:,1:nCellsSolve))
-
- !check monotonicity
- if(new_min < cur_min-eps) then
- write(*,*) 'Minimum out of bounds on tracer ', iTracer, cur_min, new_min
- end if
-
- if(new_max > cur_max+eps) then
- write(*,*) 'Maximum out of bounds on tracer ', iTracer, cur_max, new_max
- end if
- end if
- end do ! iTracer loop
-
- deallocate(tracer_new)
- deallocate(tracer_cur)
- deallocate(upwind_tendency)
- deallocate(inv_h_new)
- deallocate(tracer_max)
- deallocate(tracer_min)
- deallocate(flux_incoming)
- deallocate(flux_outgoing)
- deallocate(high_order_horiz_flux)
- deallocate(high_order_vert_flux)
- end subroutine mpas_tracer_advection_mono_tend!}}}
-
-end module mpas_tracer_advection_mono
Deleted: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -1,59 +0,0 @@
-module mpas_tracer_advection_std
-
- use mpas_kind_types
- use mpas_grid_types
- use mpas_configure
- use mpas_dmpar
- use mpas_timer
-
- 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(tracers, uh, w, verticalCellSize, grid, tend)!{{{
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! Input: s - current model state
- ! grid - grid metadata
- !
- ! Output: tend - computed scalar tendencies
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- 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) :: verticalCellSize !< Input: Distance between vertical interfaces of a cell
- type (mesh_type), intent(in) :: grid
-
- call mpas_timer_start("tracer-hadv", .false.)
- call mpas_tracer_advection_std_hadv_tend(tracers, uh, grid, tend)
- call mpas_timer_stop("tracer-hadv")
- call mpas_timer_start("tracer-vadv", .false.)
- call mpas_tracer_advection_std_vadv_tend(tracers, w, verticalCellSize, grid, tend)
- call mpas_timer_stop("tracer-vadv")
-
- end subroutine mpas_tracer_advection_std_tend!}}}
-
- subroutine mpas_tracer_advection_std_init(err)!{{{
- integer, intent(inout) :: err
-
- integer :: err_tmp
-
- err = 0
-
- 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)
-
- end subroutine mpas_tracer_advection_std_init!}}}
-
-end module mpas_tracer_advection_std
Deleted: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_hadv.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_hadv.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_hadv.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -1,105 +0,0 @@
-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%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
- tracer_weight = adv_coefs(i,iEdge) + coef_3rd_order*sign(1.0, 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
Deleted: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -1,65 +0,0 @@
-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, verticalCellSize, grid, tend)!{{{
-
- real (kind=RKIND), dimension(:,:,:), intent(inout) :: tend
- real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
- real (kind=RKIND), dimension(:,:), intent(in) :: w, verticalCellSize
- type (mesh_type), intent(in) :: grid
- real (kind=RKIND) :: dt
-
- if(order2) then
- call mpas_tracer_advection_std_vadv2_tend(tracers, w, verticalCellSize, grid, tend)
- else if(order3) then
- call mpas_tracer_advection_std_vadv3_tend(tracers, w, verticalCellSize, grid, tend)
- else if(order4) then
- call mpas_tracer_advection_std_vadv4_tend(tracers, w, verticalCellSize, 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
-
Deleted: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv2.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv2.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv2.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -1,80 +0,0 @@
-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, verticalCellSize, 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, verticalCellSize
- 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(:,:), allocatable :: vert_flux
- real (kind=RKIND) :: weightK, weightKm1
- integer :: nVertLevels, num_tracers
- integer, dimension(:), pointer :: maxLevelCell
-
- integer, parameter :: hadv_opt = 2
-
- nVertLevels = grid % nVertLevels
- num_tracers = size(tracers, dim=1)
- maxLevelCell => grid % maxLevelCell % array
-
- allocate(vert_flux(num_tracers, nVertLevels+1))
-
- !
- ! vertical flux divergence
- !
-
- ! zero fluxes at top and bottom
-
- vert_flux(:,1) = 0.
-
- do iCell=1,grid % nCellsSolve
- do k = 2, maxLevelCell(iCell)
- do iTracer=1,num_tracers
- weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*tracers(iTracer,k-1,iCell))
- end do
- end do
-
- vert_flux(:,maxLevelCell(iCell)+1) = 0
-
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + ( vert_flux(iTracer, k+1) - vert_flux(iTracer, k))
- end do
- end do
- end do
-
- deallocate(vert_flux)
-
- end subroutine mpas_tracer_advection_std_vadv2_tend!}}}
-
-end module mpas_tracer_advection_std_vadv2
Deleted: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv3.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv3.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv3.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -1,92 +0,0 @@
-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, verticalCellSize, 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, verticalCellSize
- type (mesh_type), intent(in) :: grid
- real (kind=RKIND) :: dt
-
- integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
-
- real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
- real (kind=RKIND) :: weightK, weightKm1
- integer :: nVertLevels, num_tracers
- integer, dimension(:), pointer :: maxLevelCell
-
- integer, parameter :: hadv_opt = 2
-
- coef_3rd_order = config_coef_3rd_order
-
- nVertLevels = grid % nVertLevels
- num_tracers = size(tracers, dim=1)
- maxLevelCell => grid % maxLevelCell % array
-
- allocate(vert_flux(num_tracers, nVertLevels+1))
-
- vert_flux(:,1) = 0.
-
- do iCell=1,grid % nCellsSolve
-
- k = 2
- do iTracer=1,num_tracers
- weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*tracers(iTracer,k-1,iCell))
- enddo
-
- do k=3,maxLevelCell(iCell)-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 = maxLevelCell(iCell)
-
- do iTracer=1,num_tracers
- weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*tracers(iTracer,k-1,iCell))
- enddo
-
- vert_Flux(:, maxLevelCell(iCell)+1) = 0.0
-
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (vert_flux(iTracer, k+1) - vert_flux(iTracer, k))
- end do
- end do
- end do
-
- deallocate(vert_flux)
-
- end subroutine mpas_tracer_advection_std_vadv3_tend!}}}
-
-end module mpas_tracer_advection_std_vadv3
Deleted: branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv4.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv4.F 2012-03-09 20:16:28 UTC (rev 1613)
+++ branches/ocean_projects/monotonic_advection/src/operators/mpas_tracer_advection_std_vadv4.F 2012-03-09 20:42:24 UTC (rev 1614)
@@ -1,89 +0,0 @@
-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, verticalCellSize, 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, verticalCellSize
- type (mesh_type), intent(in) :: grid
- real (kind=RKIND) :: dt
-
- integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
-
- real (kind=RKIND), dimension(:,:), allocatable :: vert_flux
- real (kind=RKIND) :: weightK, weightKm1
- integer :: nVertLevels, num_tracers
- integer, dimension(:), pointer :: maxLevelCell
-
- nVertLevels = grid % nVertLevels
- num_tracers = size(tracers, dim=1)
- maxLevelCell => grid % maxLevelCell % array
-
- allocate(vert_flux(num_tracers, nVertLevels+1))
-
- ! vertical flux divergence
- !
-
- ! zero fluxes at top and bottom
-
- vert_flux(:,1) = 0.
-
- do iCell=1,grid % nCellsSolve
-
- k = 2
- do iTracer=1,num_tracers
- weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*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 = maxLevelCell(iCell)
- do iTracer=1,num_tracers
- weightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- weightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k-1, iCell) + verticalCellSize(k, iCell))
- vert_flux(iTracer,k) = w(k,iCell)*(weightK*tracers(iTracer,k,iCell)+weightKm1*tracers(iTracer,k-1,iCell))
- enddo
-
- vert_flux(:,maxLevelCell(iCell)+1) = 0.0
-
- do k=1,maxLevelCell(iCell)
- do iTracer=1,num_tracers
- tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + (vert_flux(iTracer, k+1) - vert_flux(iTracer, k))
- 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>