<p><b>dwj07@fsu.edu</b> 2012-07-20 19:57:44 -0600 (Fri, 20 Jul 2012)</p><p><br>
-- BRANCH COMMIT --<br>
<br>
Unifing the implementation of implicit vertical mixing across time integrators.<br>
Adding halo updates of u and tracers to both, right after implicit vertical mixing is called to allow correct diagnostic<br>
values at the end of a time step.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/rich_vert_mix_bug/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- branches/ocean_projects/rich_vert_mix_bug/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-07-19 21:53:22 UTC (rev 2038)
+++ branches/ocean_projects/rich_vert_mix_bug/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-07-21 01:57:44 UTC (rev 2039)
@@ -85,11 +85,6 @@
real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
- integer :: nCells
- real (kind=RKIND), dimension(:,:), pointer :: u, h, h_edge, vertViscTopOfEdge, vertDiffTopOfCell, ke_edge
- real (kind=RKIND), dimension(:,:,:), pointer :: tracers
- integer, dimension(:), pointer :: maxLevelCell
-
block => domain % blocklist
call mpas_allocate_state(block, provis, &
block % mesh % nCells, block % mesh % nEdges, block % mesh % maxEdges, block % mesh % maxEdges2, &
@@ -262,61 +257,25 @@
! A little clean up at the end: decouple new scalar fields and compute diagnostics for new state
!
call mpas_timer_start("RK4-cleaup phase")
- block => domain % blocklist
- do while (associated(block))
- u => block % state % time_levs(2) % state % u % array
- tracers => block % state % time_levs(2) % state % tracers % array
- h => block % state % time_levs(2) % state % h % array
- h_edge => block % state % time_levs(2) % state % h_edge % array
- ke_edge => block % state % time_levs(2) % state % ke_edge % array
- vertViscTopOfEdge => block % diagnostics % vertViscTopOfEdge % array
- vertDiffTopOfCell => block % diagnostics % vertDiffTopOfCell % array
- maxLevelCell => block % mesh % maxLevelCell % array
-
- nCells = block % mesh % nCells
+ if (config_implicit_vertical_mix) then
+ call mpas_timer_start("RK4-implicit vert mix")
+ block => domain % blocklist
+ do while(associated(block))
+ call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+ block => block % next
+ end do
- do iCell=1,nCells
- do k=1,maxLevelCell(iCell)
- tracers(:,k,iCell) = tracers(:,k,iCell) / h(k,iCell)
- end do
- end do
+ ! Update halo on u and tracers, which were just updated for implicit vertical mixing. If not done,
+ ! this leads to lack of volume conservation. It is required because halo updates in RK4 are only
+ ! conducted on tendencies, not on the velocity and tracer fields. So this update is required to
+ ! communicate the change due to implicit vertical mixing across the boundary.
+ call mpas_timer_start("RK4-implicit vert mix halos")
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % u)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % tracers)
+ call mpas_timer_stop("RK4-implicit vert mix halos")
- if (config_implicit_vertical_mix) then
- call mpas_timer_start("RK4-implicit vert mix")
-
- call ocn_vmix_coefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
-
- !
- ! Implicit vertical solve for momentum
- !
- call ocn_vel_vmix_tend_implicit(block % mesh, dt, ke_edge, vertViscTopOfEdge, h, h_edge, u, err)
-
- ! mrp 110718 filter btr mode out of u
- if (config_rk_filter_btr_mode) then
- call ocn_filter_btr_mode_u(block % state % time_levs(2) % state, block % mesh)
- !block % tend % h % array(:,:) = 0.0 ! I should not need this
- endif
-
- !
- ! Implicit vertical solve for tracers
- !
-
- call ocn_tracer_vmix_tend_implicit(block % mesh, dt, vertDiffTopOfCell, h, tracers, err)
- call mpas_timer_stop("RK4-implicit vert mix")
- end if
-
- block => block % next
- end do
-
- ! Update halo on u and tracers, which were just updated for implicit vertical mixing. If not done,
- ! this leads to lack of volume conservation. It is required because halo updates in RK4 are only
- ! conducted on tendencies, not on the velocity and tracer fields. So this update is required to
- ! communicate the change due to implicit vertical mixing across the boundary.
-
- if (config_implicit_vertical_mix) then
- call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % u)
- call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % tracers)
+ call mpas_timer_stop("RK4-implicit vert mix")
end if
block => domain % blocklist
Modified: branches/ocean_projects/rich_vert_mix_bug/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- branches/ocean_projects/rich_vert_mix_bug/src/core_ocean/mpas_ocn_time_integration_split.F 2012-07-19 21:53:22 UTC (rev 2038)
+++ branches/ocean_projects/rich_vert_mix_bug/src/core_ocean/mpas_ocn_time_integration_split.F 2012-07-21 01:57:44 UTC (rev 2039)
@@ -825,37 +825,29 @@
! END large iteration loop
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- block => domain % blocklist
- do while (associated(block))
+ if (config_implicit_vertical_mix) then
+ call mpas_timer_start("se implicit vert mix")
+ block => domain % blocklist
+ do while(associated(block))
+ call ocn_vmix_implicit(dt, block % mesh, block % diagnostics, block % state % time_levs(2) % state, err)
+ block => block % next
+ end do
+ ! Update halo on u and tracers, which were just updated for implicit vertical mixing. If not done,
+ ! this leads to lack of volume conservation. It is required because halo updates in RK4 are only
+ ! conducted on tendencies, not on the velocity and tracer fields. So this update is required to
+ ! communicate the change due to implicit vertical mixing across the boundary.
+ call mpas_timer_start("se implicit vert mix halos")
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % u)
+ call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % tracers)
+ call mpas_timer_stop("se implicit vert mix halos")
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- !
- ! Implicit vertical mixing, done after timestep is complete
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ call mpas_timer_stop("se implicit vert mix")
+ end if
- u => block % state % time_levs(2) % state % u % array
- tracers => block % state % time_levs(2) % state % tracers % array
- h => block % state % time_levs(2) % state % h % array
- h_edge => block % state % time_levs(2) % state % h_edge % array
- ke_edge => block % state % time_levs(2) % state % ke_edge % array
- num_tracers = block % state % time_levs(2) % state % num_tracers
- vertViscTopOfEdge => block % diagnostics % vertViscTopOfEdge % array
- vertDiffTopOfCell => block % diagnostics % vertDiffTopOfCell % array
- maxLevelCell => block % mesh % maxLevelCell % array
- maxLevelEdgeTop => block % mesh % maxLevelEdgeTop % array
+ block => domain % blocklist
+ do while (associated(block))
- if (config_implicit_vertical_mix) then
- call ocn_vmix_coefs(block % mesh, block % state % time_levs(2) % state, block % diagnostics, err)
-
- ! Implicit vertical solve for momentum
- call ocn_vel_vmix_tend_implicit(block % mesh, dt, ke_edge, vertvisctopofedge, h, h_edge, u, err)
-
- ! Implicit vertical solve for tracers
- call ocn_tracer_vmix_tend_implicit(block % mesh, dt, vertdifftopofcell, h, tracers, err)
- end if
-
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
@@ -870,12 +862,18 @@
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
+ ! Compute velocity transport, used in advection terms of h and tracer tendency
+ block % state % time_levs(2) % state % uTransport % array(:,:) &
+ = block % state % time_levs(2) % state % u % array(:,:) &
+ + block % state % time_levs(2) % state % uBolusGM % array(:,:)
+
call mpas_reconstruct(block % mesh, block % state % time_levs(2) % state % u % array, &
- block % state % time_levs(2) % state % uReconstructX % array, &
- block % state % time_levs(2) % state % uReconstructY % array, &
- block % state % time_levs(2) % state % uReconstructZ % array, &
- block % state % time_levs(2) % state % uReconstructZonal % array, &
- block % state % time_levs(2) % state % uReconstructMeridional % array)
+ block % state % time_levs(2) % state % uReconstructX % array, &
+ block % state % time_levs(2) % state % uReconstructY % array, &
+ block % state % time_levs(2) % state % uReconstructZ % array, &
+ block % state % time_levs(2) % state % uReconstructZonal % array, &
+ block % state % time_levs(2) % state % uReconstructMeridional % array &
+ )
!TDR
call mpas_reconstruct(block % mesh, block % mesh % u_src % array, &
@@ -887,22 +885,13 @@
)
!TDR
-
call ocn_time_average_accumulate(block % state % time_levs(2) % state, block % state % time_levs(1) % state)
-
block => block % next
end do
+
call mpas_timer_stop("se timestep", timer_main)
- ! A halo update is required after implicit vertical mixing, because
- ! IVM changes the values of u and tracer, while the halo updates above
- ! are on the tendency.
- if (config_implicit_vertical_mix) then
- call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % u)
- call mpas_dmpar_exch_halo_field(domain % blocklist % state % time_levs(2) % state % tracers)
- end if
-
end subroutine ocn_time_integrator_split!}}}
end module ocn_time_integration_split
Modified: branches/ocean_projects/rich_vert_mix_bug/src/core_ocean/mpas_ocn_vmix.F
===================================================================
--- branches/ocean_projects/rich_vert_mix_bug/src/core_ocean/mpas_ocn_vmix.F 2012-07-19 21:53:22 UTC (rev 2038)
+++ branches/ocean_projects/rich_vert_mix_bug/src/core_ocean/mpas_ocn_vmix.F 2012-07-21 01:57:44 UTC (rev 2039)
@@ -576,6 +576,61 @@
!***********************************************************************
!
+! routine ocn_vmix_implicit
+!
+!> \brief Driver for implicit vertical mixing
+!> \author Doug Jacobsen
+!> \date 19 September 2011
+!> \version SVN:$Id$
+!> \details
+!> This routine is a driver for handling implicit vertical mixing
+!> of both momentum and tracers for a block. It's intended to reduce
+!> redundant code.
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_vmix_implicit(dt, grid, diagnostics, state, err)!{{{
+ real (kind=RKIND), intent(in) :: dt
+ type (mesh_type), intent(in) :: grid
+ type (diagnostic_type), intent(inout) :: diagnostics
+ type (state_type), intent(inout) :: state
+ integer, intent(out) :: err
+
+ integer :: nCells
+ real (kind=RKIND), dimension(:,:), pointer :: u, h, h_edge, vertViscTopOfEdge, vertDiffTopOfCell, ke_edge
+ real (kind=RKIND), dimension(:,:,:), pointer :: tracers
+ integer, dimension(:), pointer :: maxLevelCell
+
+ err = 0
+
+ u => state % u % array
+ tracers => state % tracers % array
+ h => state % h % array
+ h_edge => state % h_edge % array
+ ke_edge => state % ke_edge % array
+ vertViscTopOfEdge => diagnostics % vertViscTopOfEdge % array
+ vertDiffTopOfCell => diagnostics % vertDiffTopOfCell % array
+ maxLevelCell => mesh % maxLevelCell % array
+
+ nCells = mesh % nCells
+
+ call ocn_vmix_coefs(mesh, state, diagnostics, err)
+
+ !
+ ! Implicit vertical solve for momentum
+ !
+ call ocn_vel_vmix_tend_implicit(mesh, dt, ke_edge, vertViscTopOfEdge, h, h_edge, u, err)
+
+ !
+ ! Implicit vertical solve for tracers
+ !
+
+ call ocn_tracer_vmix_tend_implicit(mesh, dt, vertDiffTopOfCell, h, tracers, err)
+
+ end subroutine ocn_vmix_implicit!}}}
+
+!***********************************************************************
+!
! routine ocn_vmix_init
!
!> \brief Initializes ocean vertical mixing quantities
</font>
</pre>