<p><b>mhecht@lanl.gov</b> 2010-04-22 16:29:47 -0600 (Thu, 22 Apr 2010)</p><p>Port of O(3) advection appears to be correct now (still in development<br>
context in which it is applied to some but not all of the tracers).<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/port_adv_mwh/src/core_sw/Makefile
===================================================================
--- branches/ocean_projects/port_adv_mwh/src/core_sw/Makefile 2010-04-22 22:23:02 UTC (rev 203)
+++ branches/ocean_projects/port_adv_mwh/src/core_sw/Makefile 2010-04-22 22:29:47 UTC (rev 204)
@@ -2,6 +2,7 @@
OBJS = module_test_cases.o \
module_time_integration.o \
+ module_advection.o \
mpas_interface.o
all: core_sw
@@ -13,8 +14,10 @@
module_time_integration.o:
-mpas_interface.o: module_test_cases.o module_time_integration.o
+module_advection.o:
+mpas_interface.o: module_advection.o module_test_cases.o module_time_integration.o
+
clean:
$(RM) *.o *.mod libdycore.a
Modified: branches/ocean_projects/port_adv_mwh/src/core_sw/Registry
===================================================================
--- branches/ocean_projects/port_adv_mwh/src/core_sw/Registry 2010-04-22 22:23:02 UTC (rev 203)
+++ branches/ocean_projects/port_adv_mwh/src/core_sw/Registry 2010-04-22 22:29:47 UTC (rev 204)
@@ -7,6 +7,9 @@
namelist integer sw_model config_ntimesteps 7500
namelist integer sw_model config_output_interval 500
namelist real sw_model config_visc 0.0
+namelist integer sw_model config_tracer_adv_order 2
+namelist logical sw_model config_positive_definite false
+namelist logical sw_model config_monotonic 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
@@ -25,6 +28,8 @@
dim TWO 2
dim R3 3
dim vertexDegree vertexDegree
+dim FIFTEEN 15
+dim TWENTYONE 21
dim nVertLevels nVertLevels
dim nTracers nTracers
@@ -77,6 +82,10 @@
var real fVertex ( nVertices ) iro fVertex - -
var real h_s ( nCells ) iro h_s - -
+# Space needed for advection
+var real deriv_two ( FIFTEEN TWO nEdges ) o deriv_two - -
+var integer advCells ( TWENTYONE nCells ) - advCells - -
+
# Arrays required for reconstruction of velocity field
var real coeffs_reconstruct ( R3 maxEdges nCells ) - coeffs_reconstruct - -
Modified: branches/ocean_projects/port_adv_mwh/src/core_sw/module_time_integration.F
===================================================================
--- branches/ocean_projects/port_adv_mwh/src/core_sw/module_time_integration.F 2010-04-22 22:23:02 UTC (rev 203)
+++ branches/ocean_projects/port_adv_mwh/src/core_sw/module_time_integration.F 2010-04-22 22:29:47 UTC (rev 204)
@@ -488,24 +488,32 @@
type (grid_state), intent(in) :: s
type (grid_meta), intent(in) :: grid
- integer :: iCell, iEdge, k, iTracer, cell1, cell2
+ integer :: i, iCell, iEdge, k, iTracer, cell1, cell2
real (kind=RKIND) :: flux, tracer_edge
+ real (kind=RKIND) :: d2fdx2_cell1, d2fdx2_cell2
- !! adding pointers here
- real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell ! for areaCell, dvEdge
+ !! --- adding pointers here ---
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, areaCell
! can delete "_new"?
- real (kind=RKIND), dimension(:,:,:), pointer :: tracer_new, tracer_tend ! for tracer_tend
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracer_new, tracer_tend
integer, dimension(:,:), pointer :: cellsOnEdge
+ real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND) :: coef_3rd_order
+ dcEdge => grid % dcEdge % array
+ coef_3rd_order = 0.
+ if (config_tracer_adv_order == 3) coef_3rd_order = 1.0
+ if (config_tracer_adv_order == 3 .and. config_monotonic) coef_3rd_order = 0.25
+
+ deriv_two => grid % deriv_two % array
dvEdge => grid % dvEdge % array
tracer_new => s % tracers % array
cellsOnEdge => grid % cellsOnEdge % array
areaCell => grid % areaCell % array
tracer_tend => tend % tracers % array
- !! --------------------
+ !! ----------------------------
-
-!! okay to comment this out, relying on tend's being zeroed in _le2 routine?
+!!$! commenting this out, relying on tend's being zeroed in _le2 routine:
!!$ tracer_tend(:,:,:) = 0.0
!!$! this was the original code:
@@ -532,45 +540,71 @@
!!$ end do
!!$ end do
-!!$! here is new code:
-!!$ if (config_tracer_adv_order == 2) then
+!!$! here is new (equivalent) code:
+ if (config_tracer_adv_order == 2) then
- do iEdge=1,grid%nEdges
- cell1 = cellsOnEdge(1,iEdge)
- cell2 = cellsOnEdge(2,iEdge)
- if (cell1 > 0 .and. cell2 > 0) then
- do k=1,grid % nVertLevels
- do iTracer=3,grid % nTracers
- tracer_edge = 0.5 * (tracer_new(iTracer,k,cell1) + tracer_new(iTracer,k,cell2))
- flux = s % u % array(k,iEdge) * s % h_edge % array(k,iEdge) * dvEdge(iEdge) * tracer_edge
- tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
- tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+ do k=1,grid % nVertLevels
+ do iTracer=3,grid % nTracers
+ tracer_edge = 0.5 * (tracer_new(iTracer,k,cell1) + tracer_new(iTracer,k,cell2))
+ flux = s % u % array(k,iEdge) * s % h_edge % array(k,iEdge) * dvEdge(iEdge) * tracer_edge
+ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
+ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+ end do
end do
- end do
- end if
- end do
+ end if
+ end do
-!!$ do iEdge=1,grid%nEdges
-!!$ cell1 = cellsOnEdge(1,iEdge)
-!!$ cell2 = cellsOnEdge(2,iEdge)
-!!$ if (cell1 > 0 .and. cell2 > 0) then
-!!$ do k=1,grid % nVertLevels
-!!$ do iTracer=3,grid % nTracers
-!!$ tracer_edge = 0.5 * (tracer_new(iTracer,k,cell1) + tracer_new(iTracer,k,cell2))
-!!$ flux = s % u % array(k,iEdge) * s % h_edge % array(k,iEdge) * dvEdge(iEdge) * tracer_edge
-!!$ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
-!!$ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
-!!$ end do
-!!$ end do
-!!$ end if
-!!$ end do
+ else if (config_tracer_adv_order == 3) then
-!!$ else if (config_tracer_adv_order == 3) then
-!!$ endif
+ do iEdge=1,grid%nEdges
+ cell1 = cellsOnEdge(1,iEdge)
+ cell2 = cellsOnEdge(2,iEdge)
+ if (cell1 > 0 .and. cell2 > 0) then
+ do k=1,grid % nVertLevels
+
+ do iTracer=3,grid % nTracers
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * tracer_new(iTracer,k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * tracer_new(iTracer,k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + &
+ deriv_two(i+1,1,iEdge) * tracer_new(iTracer,k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + &
+ deriv_two(i+1,2,iEdge) * tracer_new(iTracer,k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ if (s % u % array(k,iEdge) * s % h_edge % array(k,iEdge) > 0) then
+ flux = dvEdge(iEdge) * s % u % array(k,iEdge) * s % h_edge % array(k,iEdge) * ( &
+ 0.5*(tracer_new(iTracer,k,cell1) + tracer_new(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ -(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+ else
+ flux = dvEdge(iEdge) * s % u % array(k,iEdge) * s % h_edge % array(k,iEdge) * ( &
+ 0.5*(tracer_new(iTracer,k,cell1) + tracer_new(iTracer,k,cell2)) &
+ -(dcEdge(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12. &
+ +(dcEdge(iEdge) **2) * coef_3rd_order*(d2fdx2_cell1 - d2fdx2_cell2) / 12. )
+ end if
+
+ tracer_tend(iTracer,k,cell1) = tracer_tend(iTracer,k,cell1) - flux/areaCell(cell1)
+ tracer_tend(iTracer,k,cell2) = tracer_tend(iTracer,k,cell2) + flux/areaCell(cell2)
+
+ end do
+ end do
+ end if
+ end do
+
+ endif
+
end subroutine compute_scalar_tend_gt2
-
subroutine compute_solve_diagnostics(dt, s, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields used in the tendency computations
Modified: branches/ocean_projects/port_adv_mwh/src/core_sw/mpas_interface.F
===================================================================
--- branches/ocean_projects/port_adv_mwh/src/core_sw/mpas_interface.F 2010-04-22 22:23:02 UTC (rev 203)
+++ branches/ocean_projects/port_adv_mwh/src/core_sw/mpas_interface.F 2010-04-22 22:29:47 UTC (rev 204)
@@ -16,6 +16,7 @@
subroutine mpas_init(block, mesh, dt)
use grid_types
+ use advection
use time_integration
implicit none
@@ -27,6 +28,7 @@
call compute_solve_diagnostics(dt, block % time_levs(1) % state, mesh)
call init_reconstruct(mesh)
call reconstruct(block % time_levs(1) % state, mesh)
+ call initialize_advection_rk(mesh)
end subroutine mpas_init
</font>
</pre>