<p><b>dwj07@fsu.edu</b> 2012-07-23 14:57:17 -0600 (Mon, 23 Jul 2012)</p><p><br>
-- TRUNK COMMIT --<br>
<br>
core_ocean only<br>
Adding a halo update for implicit vertical mixing.<br>
Fixing issues with richardson number vertical mixing.<br>
</p><hr noshade><pre><font color="gray">Modified: trunk/mpas/src/core_ocean/mpas_ocn_equation_of_state.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_equation_of_state.F 2012-07-23 20:32:07 UTC (rev 2044)
+++ trunk/mpas/src/core_ocean/mpas_ocn_equation_of_state.F 2012-07-23 20:57:17 UTC (rev 2045)
@@ -74,10 +74,13 @@
! Input: grid - grid metadata
! s - state: tracers
! k_displaced
- ! If k_displaced<=0, state % rho is returned with no displaced
- ! If k_displaced>0,the state % rhoDisplaced is returned, and is for
+ !
+ ! If k_displaced==0, state % rho is returned with no displacement
+ !
+ ! If k_displaced~=0,the state % rhoDisplaced is returned, and is for
! a parcel adiabatically displaced from its original level to level
- ! k_displaced. This does not effect the linear EOS.
+ ! k_displaced. When using the linear EOS, state % rhoDisplaced is
+ ! still filled, but depth (i.e. pressure) does not modify the output.
!
! Output: s - state: computed density
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -103,19 +106,19 @@
indexT = s % index_temperature
indexS = s % index_salinity
- if (linearEos) then
+ ! Choose to fill the array rho or rhoDisplaced
+ if (k_displaced == 0) then
rho => s % rho % array
+ else
+ rho => s % rhoDisplaced % array
+ endif
+ if (linearEos) then
+
call ocn_equation_of_state_linear_rho(grid, indexT, indexS, tracers, rho, err)
elseif (jmEos) then
- if(k_displaced == 0) then
- rho => s % rho % array
- else
- rho => s % rhoDisplaced % array
- endif
-
call ocn_equation_of_state_jm_rho(grid, k_displaced, displacement_type, indexT, indexS, tracers, rho, err)
endif
Modified: trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2012-07-23 20:32:07 UTC (rev 2044)
+++ trunk/mpas/src/core_ocean/mpas_ocn_mpas_core.F 2012-07-23 20:57:17 UTC (rev 2045)
@@ -259,7 +259,7 @@
call ocn_diagnostic_solve(dt, block % state % time_levs(1) % state, mesh)
call mpas_timer_stop("diagnostic solve", initDiagSolveTimer)
- ! Compute velocity transport, used in advection terms of h and tracer tendancy
+ ! Compute velocity transport, used in advection terms of h and tracer tendency
block % state % time_levs(1) % state % uTransport % array(:,:) &
= block % state % time_levs(1) % state % u % array(:,:) &
+ block % state % time_levs(1) % state % uBolusGM % array(:,:)
Modified: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-07-23 20:32:07 UTC (rev 2044)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_rk4.F 2012-07-23 20:57:17 UTC (rev 2045)
@@ -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, &
@@ -217,7 +212,7 @@
call ocn_diagnostic_solve(dt, provis, block % mesh)
- ! Compute velocity transport, used in advection terms of h and tracer tendancy
+ ! Compute velocity transport, used in advection terms of h and tracer tendency
provis % uTransport % array(:,:) &
= provis % u % array(:,:) &
+ provis % uBolusGM % array(:,:)
@@ -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 weres 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
@@ -336,7 +295,7 @@
call ocn_diagnostic_solve(dt, block % state % time_levs(2) % state, block % mesh)
- ! Compute velocity transport, used in advection terms of h and tracer tendancy
+ ! 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(:,:)
Modified: trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2012-07-23 20:32:07 UTC (rev 2044)
+++ trunk/mpas/src/core_ocean/mpas_ocn_time_integration_split.F 2012-07-23 20:57:17 UTC (rev 2045)
@@ -293,7 +293,7 @@
! uTranport = uBcl + uBolus
! This is u used in advective terms for h and tracers
- ! in tendancy calls in stage 3.
+ ! in tendency calls in stage 3.
block % state % time_levs(2) % state % uTransport % array(k,iEdge) &
= block % mesh % edgeMask % array(k,iEdge) &
*( block % state % time_levs(2) % state % uBcl % array(k,iEdge) &
@@ -639,7 +639,7 @@
! uTranport = uBtr + uBcl + uBolus + uCorrection
! This is u used in advective terms for h and tracers
- ! in tendancy calls in stage 3.
+ ! in tendency calls in stage 3.
block % state % time_levs(2) % state % uTransport % array(k,iEdge) &
= block % mesh % edgeMask % array(k,iEdge) &
*( block % state % time_levs(2) % state % uBtr % array( iEdge) &
@@ -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 stage 3 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,12 +885,11 @@
)
!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)
end subroutine ocn_time_integrator_split!}}}
Modified: trunk/mpas/src/core_ocean/mpas_ocn_vmix.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_vmix.F 2012-07-23 20:32:07 UTC (rev 2044)
+++ trunk/mpas/src/core_ocean/mpas_ocn_vmix.F 2012-07-23 20:57:17 UTC (rev 2045)
@@ -47,7 +47,8 @@
ocn_tracer_vmix_tend_explicit, &
ocn_vel_vmix_tend_implicit, &
ocn_tracer_vmix_tend_implicit, &
- ocn_vmix_init
+ ocn_vmix_init, &
+ ocn_vmix_implicit
!--------------------------------------------------------------------
!
@@ -576,6 +577,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 (diagnostics_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 => grid % maxLevelCell % array
+
+ nCells = grid % nCells
+
+ call ocn_vmix_coefs(grid, state, diagnostics, err)
+
+ !
+ ! Implicit vertical solve for momentum
+ !
+ call ocn_vel_vmix_tend_implicit(grid, dt, ke_edge, vertViscTopOfEdge, h, h_edge, u, err)
+
+ !
+ ! Implicit vertical solve for tracers
+ !
+
+ call ocn_tracer_vmix_tend_implicit(grid, dt, vertDiffTopOfCell, h, tracers, err)
+
+ end subroutine ocn_vmix_implicit!}}}
+
+!***********************************************************************
+!
! routine ocn_vmix_init
!
!> \brief Initializes ocean vertical mixing quantities
Modified: trunk/mpas/src/core_ocean/mpas_ocn_vmix_coefs_rich.F
===================================================================
--- trunk/mpas/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-07-23 20:32:07 UTC (rev 2044)
+++ trunk/mpas/src/core_ocean/mpas_ocn_vmix_coefs_rich.F 2012-07-23 20:57:17 UTC (rev 2045)
@@ -472,7 +472,7 @@
drhoTopOfCell = 0.0
do iCell=1,nCells
do k=2,maxLevelCell(iCell)
- drhoTopOfCell(k,iCell) = rho(k-1,iCell) - rhoDisplaced(k-1,iCell)
+ drhoTopOfCell(k,iCell) = rhoDisplaced(k-1,iCell) - rhoDisplaced(k,iCell)
end do
end do
</font>
</pre>