<p><b>dwj07@fsu.edu</b> 2012-03-14 15:48:15 -0600 (Wed, 14 Mar 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Adding capabilities for 2nd, 3rd, and 4th order standard advection.<br>
<br>
Also, adding flags for prescribed velocity and thickness.<br>
<br>
Cleaning up some of the arrays required for advection.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/Registry 2012-03-14 21:21:26 UTC (rev 1639)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/Registry 2012-03-14 21:48:15 UTC (rev 1640)
@@ -11,6 +11,8 @@
namelist character sw_model config_run_duration none
namelist integer sw_model config_stats_interval 100
namelist logical sw_model config_initial_stats false
+namelist logical sw_model config_prescribe_velocity false
+namelist logical sw_model config_prescribe_thickness false
namelist character io config_input_name grid.nc
namelist character io config_output_name output.nc
namelist character io config_restart_name restart.nc
@@ -69,7 +71,7 @@
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 real advection config_coef_3rd_order 0.25
namelist logical advection config_monotonic false
namelist logical advection config_check_monotonicity false
namelist logical restore config_restoreTS false
@@ -83,6 +85,7 @@
dim nEdges nEdges
dim maxEdges maxEdges
dim maxEdges2 maxEdges2
+dim nAdvectionCells maxEdges2+0
dim nVertices nVertices
dim TWO 2
dim R3 3
@@ -150,14 +153,13 @@
var persistent real h_s ( nCells ) 0 iro h_s mesh - -
% Space needed for advection
-var persistent real deriv_two ( FIFTEEN TWO nEdges ) 0 - deriv_two mesh - -
-var persistent integer advCells ( TWENTYONE nCells ) 0 - advCells mesh - -
+var persistent real deriv_two ( maxEdges2 TWO nEdges ) 0 - deriv_two mesh - -
% Added for monotonic advection scheme
-var persistent real adv_coefs ( FIFTEEN nEdges ) 0 - adv_coefs mesh - -
-var persistent real adv_coefs_2nd ( FIFTEEN nEdges ) 0 - adv_coefs_2nd mesh - -
-var persistent real adv_coefs_3rd ( FIFTEEN nEdges ) 0 - adv_coefs_3rd mesh - -
-var persistent integer advCellsForEdge ( FIFTEEN nEdges ) 0 - advCellsForEdge mesh - -
+var persistent real adv_coefs ( nAdvectionCells nEdges ) 0 - adv_coefs mesh - -
+var persistent real adv_coefs_2nd ( nAdvectionCells nEdges ) 0 - adv_coefs_2nd mesh - -
+var persistent real adv_coefs_3rd ( nAdvectionCells nEdges ) 0 - adv_coefs_3rd mesh - -
+var persistent integer advCellsForEdge ( nAdvectionCells nEdges ) 0 - advCellsForEdge mesh - -
var persistent integer nAdvCellsForEdge ( nEdges ) 0 - nAdvCellsForEdge mesh - -
var persistent integer highOrderAdvectionMask ( nVertLevels nEdges ) 0 - highOrderAdvectionMask mesh - -
var persistent integer lowOrderAdvectionMask ( nVertLevels nEdges ) 0 - lowOrderAdvectionMask mesh - -
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_advection.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_advection.F 2012-03-14 21:21:26 UTC (rev 1639)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_advection.F 2012-03-14 21:48:15 UTC (rev 1640)
@@ -27,7 +27,7 @@
integer, intent(out) :: err
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
- integer, dimension(:,:), pointer :: advCells
+ integer, dimension(:), pointer :: advCells
! local variables
@@ -76,7 +76,8 @@
pii = 2.*asin(1.0)
- advCells => grid % advCells % array
+! advCells => grid % advCells % array
+ allocate(advCells(grid % maxEdges2))
deriv_two => grid % deriv_two % array
deriv_two(:,:,:) = 0.
@@ -104,7 +105,7 @@
end do
end if
- advCells(1,iCell) = n
+ advCells(1) = n
! check to see if we are reaching outside the halo
@@ -121,10 +122,10 @@
if ( grid % on_a_sphere ) then
do i=1,n
- advCells(i+1,iCell) = cell_list(i)
- xc(i) = grid % xCell % array(advCells(i+1,iCell))/a
- yc(i) = grid % yCell % array(advCells(i+1,iCell))/a
- zc(i) = grid % zCell % array(advCells(i+1,iCell))/a
+ advCells(i+1) = cell_list(i)
+ xc(i) = grid % xCell % array(advCells(i+1))/a
+ yc(i) = grid % yCell % array(advCells(i+1))/a
+ zc(i) = grid % zCell % array(advCells(i+1))/a
end do
theta_abs(iCell) = pii/2. - sphere_angle( xc(1), yc(1), zc(1), &
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-14 21:21:26 UTC (rev 1639)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_mpas_core.F 2012-03-14 21:48:15 UTC (rev 1640)
@@ -731,8 +731,12 @@
if (boundaryEdge(k,iEdge).eq.1) then
boundaryCell(k,cellsOnEdge(1,iEdge)) = 1
boundaryCell(k,cellsOnEdge(2,iEdge)) = 1
- cellMask(k,cellsOnEdge(1,iEdge)) = 0
- cellMask(k,cellsOnEdge(2,iEdge)) = 0
+ if(maxLevelCell(cellsOnEdge(1,iEdge)) > k) then
+ cellMask(k, cellsOnEdge(1,iEdge)) = 0
+ end if
+ if(maxLevelCell(cellsOnEdge(2,iEdge)) > k) then
+ cellMask(k, cellsOnEdge(2,iEdge)) = 0
+ end if
boundaryVertex(k,verticesOnEdge(1,iEdge)) = 1
boundaryVertex(k,verticesOnEdge(2,iEdge)) = 1
vertexMask(k,verticesOnEdge(1,iEdge)) = 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-14 21:21:26 UTC (rev 1639)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-03-14 21:48:15 UTC (rev 1640)
@@ -231,6 +231,14 @@
provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
+ if (config_prescribe_velocity) then
+ provis % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ end if
+
+ if (config_prescribe_thickness) then
+ provis % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ end if
+
call ocn_diagnostic_solve(dt, provis, block % mesh)
block => block % next
@@ -325,6 +333,14 @@
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
+ if (config_prescribe_velocity) then
+ block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ end if
+
+ if (config_prescribe_thickness) then
+ block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ end if
+
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
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-14 21:21:26 UTC (rev 1639)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_time_integration_split.F 2012-03-14 21:48:15 UTC (rev 1640)
@@ -900,6 +900,15 @@
if (config_test_case == 1) then ! For case 1, wind field should be fixed
block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
end if
+
+ if (config_prescribe_velocity) then
+ block % state % time_levs(2) % state % u % array(:,:) = block % state % time_levs(1) % state % u % array(:,:)
+ end if
+
+ if (config_prescribe_thickness) then
+ block % state % time_levs(2) % state % h % array(:,:) = block % state % time_levs(1) % state % h % array(:,:)
+ end if
+
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
block % state % time_levs(2) % state % uReconstructX % array, &
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection.F 2012-03-14 21:21:26 UTC (rev 1639)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection.F 2012-03-14 21:48:15 UTC (rev 1640)
@@ -85,6 +85,9 @@
highOrderAdvectionMask = 0
lowOrderAdvectionMask = 0
+ if(config_horiz_tracer_adv_order == 2) then
+
+ end if
do iEdge = 1, grid % nEdges
nAdvCellsForEdge(iEdge) = 0
@@ -231,6 +234,12 @@
deallocate(cell_list)
deallocate(ordered_cell_list)
+ ! If 2nd order advection, set masks appropriately.
+ if(config_horiz_tracer_adv_order == 2) then
+ lowOrderAdvectionMask = 1
+ highOrderAdvectionMask = 0
+ end if
+
if (maxval(highOrderAdvectionMask+lowOrderAdvectionMask) > 1) then
write(*,*) "Masks don't sum to 1."
err = 1
@@ -292,6 +301,7 @@
err = 0
call mpas_ocn_tracer_advection_std_init(err_tmp)
+ call mpas_ocn_tracer_advection_mono_init(err_tmp)
err = ior(err, err_tmp)
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_mono.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-03-14 21:21:26 UTC (rev 1639)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_mono.F 2012-03-14 21:48:15 UTC (rev 1640)
@@ -23,8 +23,11 @@
private
save
- public :: mpas_ocn_tracer_advection_mono_tend
+ real (kind=RKIND) :: coef_3rd_order
+ public :: mpas_ocn_tracer_advection_mono_tend, &
+ mpas_ocn_tracer_advection_mono_init
+
contains
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
@@ -63,7 +66,7 @@
integer, dimension(:), pointer :: nAdvCellsForEdge, maxLevelCell, maxLevelEdgeTop, nEdgesOnCell
integer, dimension(:,:), pointer :: cellsOnEdge, cellsOnCell, advCellsForEdge, highOrderAdvectionMask, lowOrderAdvectionMask
- real (kind=RKIND) :: coef_3rd_order, flux_upwind, tracer_min_new, tracer_max_new, tracer_upwind_new, scale_factor
+ real (kind=RKIND) :: 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
@@ -74,8 +77,6 @@
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
@@ -351,4 +352,31 @@
deallocate(high_order_vert_flux)
end subroutine mpas_ocn_tracer_advection_mono_tend!}}}
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! routine mpas_ocn_tracer_advection_mono_init
+!
+!> \brief MPAS ocean initialize monotonic tracer advection tendency with FCT
+!> \author Doug Jacobsen
+!> \date 03/09/12
+!> \version SVN:$Id:$
+!> \details
+!> This routine initializes the monotonic tracer advection tendencity using a FCT.
+!
+!-----------------------------------------------------------------------
+ subroutine mpas_ocn_tracer_advection_mono_init(err)!{{{
+ integer, intent(inout) :: err !< Input: Error Flags
+
+ integer :: err_tmp
+
+ err = 0
+
+ if ( config_horiz_tracer_adv_order == 3) then
+ coef_3rd_order = config_coef_3rd_order
+ else if(config_horiz_tracer_adv_order == 2 .or. config_horiz_tracer_adv_order == 4) then
+ coef_3rd_order = 0.0
+ end if
+
+ end subroutine mpas_ocn_tracer_advection_mono_init!}}}
+
end module mpas_ocn_tracer_advection_mono
Modified: branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_std_hadv.F
===================================================================
--- branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_std_hadv.F 2012-03-14 21:21:26 UTC (rev 1639)
+++ branches/ocean_projects/monotonic_advection/src/core_ocean/mpas_ocn_tracer_advection_std_hadv.F 2012-03-14 21:48:15 UTC (rev 1640)
@@ -57,9 +57,9 @@
real (kind=RKIND), dimension(:), pointer :: areaCell
integer, dimension(:,:), pointer :: cellsOnEdge
- integer, dimension(:,:), pointer :: advCellsForEdge
+ integer, dimension(:,:), pointer :: advCellsForEdge, highOrderAdvectionMask, lowOrderAdvectionMask
integer, dimension(:), pointer :: nAdvCellsForEdge
- real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_3rd
+ real (kind=RKIND), dimension(:,:), pointer :: adv_coefs, adv_coefs_2nd, adv_coefs_3rd
real (kind=RKIND), dimension(:,:), allocatable :: flux_arr
integer :: nVertLevels, num_tracers
@@ -71,7 +71,10 @@
nAdvCellsForEdge => grid % nAdvCellsForEdge % array
advCellsForEdge => grid % advCellsForEdge % array
adv_coefs => grid % adv_coefs % array
+ adv_coefs_2nd => grid % adv_coefs_2nd % array
adv_coefs_3rd => grid % adv_coefs_3rd % array
+ highOrderAdvectionMask => grid % highOrderAdvectionMask % array
+ lowOrderAdvectionMask => grid % lowOrderAdvectionMask % array
nVertLevels = grid % nVertLevels
num_tracers = size(tracers, dim=1)
@@ -87,7 +90,8 @@
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)
+ tracer_weight = lowOrderAdvectionMask(k, iEdge) * adv_coefs_2nd(i,iEdge) &
+ + highOrderAdvectionMask(k, iEdge) * (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
@@ -126,7 +130,11 @@
hadvOn =.true.
- coef_3rd_order = config_coef_3rd_order
+ if ( config_horiz_tracer_adv_order == 3) then
+ coef_3rd_order = config_coef_3rd_order
+ else if ( config_horiz_tracer_adv_order == 2 .or. config_horiz_tracer_adv_order == 4) then
+ coef_3rd_order = 0.0
+ end if
end subroutine mpas_ocn_tracer_advection_std_hadv_init!}}}
end module mpas_ocn_tracer_advection_std_hadv
</font>
</pre>